function [p,dr,ps,de,der]=getinv(di,p,ps,dr,iplot)

% function [p,dr,ps,de,der]=getinv(di,p,ps,dr,iplot)

% LADCP processing software version 5.0

%

% - solve linear inverse problem

% 

%  Martin Visbeck, LDEO, April-2000

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

%                    G E T I N V . M 

%                    doc: Thu Jun 17 15:36:21 2004

%                    dlm: Wed Sep  4 16:55:50 2019

%                    (c) 2004 ladcp@

%                    uE-Info: 35 55 NIL 0 0 72 0 2 8 NIL ofnI

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

% CHANGE HISTORY:

%

% Jun 17, 2004: - removed ps.botfac = 0 on lack of magdecl

% Jun 27, 2004: - added warning about removed SADCP profiles

% Jul  6, 2004: - BUG: warning did not work because p was not accessible

%		       inside lainsadcp()

% Jul 17, 2004: - added warning about down/up mismatch (GPS problems)

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

% Aug 25, 2010: - BUG: warning about low SADCP weight was produced

%		       (and corresponding velocities removed) even

%		       on sadcpfac=0, i.e. when SADCP data were only

%		       used for plotting

% Jun 29, 2011: - pass ps to geterr to allow controlling fig. 3

% Jan  6, 2012: - lesqchol removed because it does not deal with

%		  nearly singular matrices gracefully

%		- BUG: replaced all imag() by new imagnan()

% Jul 28, 2014: - modified how to specify smallfac

% Aug  9, 2014: - modified how to specify dragfac (+ve for fixed; -ve for scaled Martin's default)

% Jul 23, 2015: - commented on bug below (#!#)

%		- rounded default ps.dz (GK suggestion)

if nargin<5, iplot=0; end

%###  start by defining all processing parameter which are not set by the user

disp('GETINV:  compute best velocity profile')

ps=setdefv(ps,'dz',round(medianan(abs(diff(di.izm(:,1))))));

ps=setdefv(ps,'botfac',1); 

ps=setdefv(ps,'sadcpfac',1); 

ps=setdefv(ps,'barofac',1);

ps=setdefv(ps,'dragfac',0);

% negative values weigh Martin's original default

if (ps.dragfac < 0)

  ps.dragfac = -ps.dragfac * tanh(p.maxdepth/3000)*0.2;

end

% how much do you want to smooth the ocean profile

ps=setdefv(ps,'smoofac',0);

% how much do you want to restrict large shears on short vertical wavelengths

%  	ps.smallfac(1) is the vertical-wavelength filter cutoff

%	ps.smallfac(2) is the filter strength (larger implies more filtering)

%ps=setdefv(ps,'smallfac',[500 0.02]);		% Martin's defaults

%ps=setdefv(ps,'smallfac',[200 0.02]);		% works reasonably well for DoMORE1 data

ps=setdefv(ps,'smallfac',[99999 0]);		% default: disable

    imax = ceil(p.maxdepth/ps.smallfac(1));

    smallfac = [0 0];

    for i=1:imax

	smallfac(i,1) = i;

	smallfac(i,2) = ps.smallfac(2)/(1+abs(i-(imax/2)))*tanh(p.maxdepth/3000);

    end

% do you want up/down cast seperately then set to 1

ps=setdefv(ps,'down_up',1);

% decide how to weight data

ps=setdefv(ps,'std_weight',1);

% taper small weights weight = weight. ^weightpower

ps=setdefv(ps,'weightpower',1); 

ps=setdefv(ps,'solve',1); 

% weight bottom track data with distance of bottom

% ps=setdefv(ps,'btrk_weight_nblen',[15 5]);

ps=setdefv(ps,'btrk_weight_nblen',[0 0]);

% restrict data to one instrument only 1: up+dn, 2:dn only  3:up only

ps=setdefv(ps,'up_dn_looker',1);

% Check for magnetic deviation

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;

  p.rot=0;

  ps.sadcpfac=0;

  ps.barofac=0;

%  ps.botfac=0; %%% NOT NEEDED?!

  ps.dragfac=0;

 end

else

  warn=[' NO magnetic deviation given '];

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

  p.rot=0;

  ps.sadcpfac=0;

  ps.barofac=0;

%  ps.botfac=0; %%% NOT NEEDED?!

  ps.dragfac=0;

end

% Barotropic velocity error due to navigation error

ps=setdefv(ps,'barvelerr',2*p.nav_error/p.dt_profile);

disp([' Barotropic velocity error ',num2str(2*p.nav_error/p.dt_profile),...

      ' [m/s]'])

% Super ensemble velocity error

nmax=min(length(di.izd),7);

sw=stdnan(di.rw(di.izd(1:nmax),:)); ii=find(sw>0);

sw=medianan(sw(ii))/tan(p.beamangle*pi/180);

disp([' super ensemble velocity error ',num2str(sw),' [m/s]'])

ps=setdefv(ps,'velerr',maxnan([sw,0.02])); 

if exist('dr','var')

 if existf(dr,'uerr')

  if any(isfinite(dr.uerr))

   ps.velerr=medianan(dr.uerr);

   disp([' set velocity error to:',num2str(ps.velerr),' [m/s]'])

  end

 end

end

disp([' vertical resolution (ps.dz) is ',num2str(ps.dz),' [m]'])

[nbin,nt]=size(di.ru);

%### set up matrixes for inverse problem

%set up data arrays 

% restrict to up/down looker

if ps.up_dn_looker==2

 di.weight(di.izu,:)=nan;

 disp(' restrict inversion to down looking instrument only')

elseif ps.up_dn_looker==3

 di.weight(di.izd,:)=nan;

 disp(' restrict inversion to up looking instrument only')

 p=rmfield(p,'zbottom');

end

%  velocities are complex

d=di.ru+sqrt(-1)*di.rv;

d=reshape(d,nbin*nt,1);

% bottom track velocity

if existf(p,'zbottom');

 bvel=di.bvel(1,:)+sqrt(-1)*di.bvel(2,:);

 bvels=sqrt(di.bvels(1,:).^2+di.bvels(2,:).^2);

 if sum(ps.btrk_weight_nblen) > 0

  % normalize btweight with gaussian centered at 5 binlenght from bottom

  hbot=di.z+p.zbottom;

  hm=abs(diff(di.izm([-1:0]+end,1)))*ps.btrk_weight_nblen(1);

  hs=abs(diff(di.izm([-1:0]+end,1)))*ps.btrk_weight_nblen(2);

  whbot=exp(-(hbot-hm).^2./(hs).^2);

  bvels=bvels./whbot;

  disp([' weighted bottom track with Gaussian centered ',...

        int2str(ps.btrk_weight_nblen(1)),' bins above bottom and width of ',...

        int2str(ps.btrk_weight_nblen(2)),' bins '])

 end

 dbot=meshgrid(bvel,1:nbin);

 dbot=reshape(dbot,nt*nbin,1);

end

% bin depths

izv=reshape(-di.izm,nbin*nt,1);

% profile number

jprof=cumsum(ones(nbin,nt)')';

jprof=reshape(jprof,nbin*nt,1);

% bin number

jbin=meshgrid(1:nbin,1:nt)';

jbin=reshape(jbin,nbin*nt,1);

% data weight

if ps.std_weight~=1;

 disp(' use correlation based weights')

 % use correlation based estimator

 wm=di.weight.^ps.weightpower;

 % data with a d.weight smaller that weightmin are not used

 ps=setdefv(ps,'weightmin',0.05); 

else

 % use super ensemble based estimator

 disp([' use super ensemble std based weights normalized by ',...

       num2str(ps.velerr),' m/s '])

 wm=ps.velerr./di.ruvs+di.weight*0;

 % data with a d.weight smaller that weightmin are not used

 ps=setdefv(ps,'weightmin',0.05); 

end

wm=reshape(wm,nt*nbin,1);

% other arrays used

dtiv=di.dt;

tim=di.time_jul;

zctd=di.z;

slat=di.slat;

slon=di.slon;

if isfinite(sum(slon+slat))

 xship=(slon-slon(1))*60*1852*cos(meannan(slat)*pi/180);

 yship=(slat-slat(1))*60*1852;

 uship=diff(xship)./diff(tim)/(24*3600);

 vship=diff(yship)./diff(tim)/(24*3600);

 shipvel=uship+sqrt(-1)*vship;

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

% get a smooth ship surface velocity

 tim1=tim-tim(1);

 shipvelf=polyval(polyfit(tim1,shipvel,3),tim1);

 i=0;

 if length(shipvel)>5

  while (i<50 & std(abs(shipvel-shipvelf))>0.15) | i==0

   shipvel=mean([shipvel([3:end,end,end]);shipvel([2:end,end]);...

         shipvel;shipvel([1,1:(end-1)]);shipvel([1,1,1:(end-2)])]);

   i=i+1;

  end

  disp([' preaveraged GPS ships vel ',int2str(i),' times '])

 end

else

 % no GPS ship navigation time series, assume constant ships velocity

 uship_a=p.uship+sqrt(-1)*p.vship;

 shipvel=uship_a+di.z*0;

end

% mean CTD vertical velocity

wctd=meannan(di.rw);

if ps.dragfac>0

  shipdragvel=shipvel;

  if exist('dr','var')

  % compute the ocean flow at the depth of the CTD

   ut=interp1(dr.z,dr.u,-di.z);

   vt=interp1(dr.z,dr.v,-di.z);

   shipdragvel=shipvel-(ut+sqrt(-1)*vt);

  end

  % derive estimate for CTD velocity using ship velocity

  % and vertical velocity and wireangle

  %  0.2 of W-velocity project horizontal for 1 m/s

  ps=setdefv(ps,'drag_tilt_vel',0.5);

  tiltfac = ps.drag_tilt_vel*abs(shipdragvel).*(1-tanh(-di.z/8000));

  % preject fration of vertical velocity in ships velocity direction

  wctdf=-shipdragvel./(abs(shipdragvel)+0.001).*wctd.*tiltfac;

  % make sure that this is not faster than ship velocity

  wfac=max(abs(wctdf)./(abs(shipdragvel)+0.001),1+wctd*0);

  wctdf=wctdf./wfac;

  % average ctdvel back in time to take the inertia of wire into account

  % smooth over max of 20 minutes for depth ~ 2000m

  ps=setdefv(ps,'drag_lag_tim',20);

  ps=setdefv(ps,'drag_lag_depth',2000);

  for j=1:length(shipdragvel)

    zfac=tanh(-di.z(j)/ps.drag_lag_depth);

    % how many super ensembles need to be averaged

    nsmooth=sum(di.time_jul>(di.time_jul(j)-ps.drag_lag_tim/24/60) & di.time_jul<di.time_jul(j));

    ii=j-[0:fix(nsmooth*zfac)];

    ii=ii(find(ii>0));

    if length(ii)<2

     ii=j;

    end

    % smooth w less in time

    iw=ii(1:fix(end/10));

    if length(iw)<2

     iw=j;

    end

    % set the assumed ctd velocity as the sum between ship vel and projected w-vel

    ctdvel(j)=meannan(shipvel(ii))+meannan(wctdf(iw));

  end

else

  ctdvel=shipvel*NaN;

end

% remove last data bin  and one bin off the bottom and surface 

if existf(p,'zbottom'),

 zbottom=p.zbottom;

else

 zbottom=1e32;

end

[jmax,jbott]=max(izv);

%#!# the following line is wrong, as izv>(ps.dz*(jmax-1) is never true;

%#!# the code is supposed to remove the "last data bin" (whatever that means)

%#!# but it does not do anything

ii=find(izv<(ps.dz) | izv>(ps.dz*(jmax-1)) | izv>(zbottom-ps.dz));

wm(ii)=0;

% remove bad or weakly constained data

jm=max(jprof);

disp([' remove ',int2str(sum(wm<ps.weightmin)),' constaints below minimum weight ']);

ii=find(isnan(d) | isnan(wm) | wm<ps.weightmin);

d(ii)=[];

if exist('dbot','var')

 dbot(ii)=[];

end

izv(ii)=[];

jprof(ii)=[];

jbin(ii)=[];

wm(ii)=[];

% remove empty profiles at the end

ii=(max(jprof)+1):jm;

tim(ii)=[];

ctdvel(ii)=[];

shipvel(ii)=[];

zctd(ii)=[];

wctd(ii)=[];

slat(ii)=[];

slon(ii)=[];

dtiv(ii)=[];

if exist('bvel','var')

 bvel(ii)=[];

 bvels(ii)=[];

end

[jmax,jbott]=max(izv);

% prepare some output arrays

if existf(p,'name')

 dr.name=p.name;

end

dr.date=gregoria(median(tim));

dr.lat=p.lat;

dr.lon=p.lon;

% set up main matrices for inversion

A1=lainseta(jprof,1);

A2=lainseta(izv,ps.dz);

[nt,nz]=size(A2);

% resulting depth vector

%z =([1:nz]'-.5)*ps.dz;

z =([1:nz]')*ps.dz;

A1o=A1;

A2o=A2;

do=d;

%### add weights to data

[A2,A1,d,idoc,iupc]=lainweig(A2,A1,d,wm);

% make sure time and depth dimension are different

if size(A2,2)==size(A1,2)

 disp('change dimension of A2')

 A2=[A2,A2(:,1)*0]; 

 A2o=[A2o,A2o(:,1)*0];

 [nt,nz]=size(A2);

 z =([1:nz]')*ps.dz;

end

% save constraints

de.ocean_constraints=full(sum(abs(A2)));

de.ctd_constraints=full(sum(abs(A1)));

de.type_constraints='Velocity  ';

%### smooth ocean and CTD velocity profiles

disp(' smooth Ocean velocity profile')

[A2,A1,d]=lainsmoo(A2,A1,d,ps.smoofac);

disp(' smooth CTD velocity profile')

[A1,A2,d]=lainsmoo(A1,A2,d,ps.smoofac);

de.ocean_constraints=[de.ocean_constraints;sum(abs(A2))-(de.ocean_constraints)];

de.ctd_constraints=[de.ctd_constraints;sum(abs(A1))-(de.ctd_constraints)];

de.type_constraints=[de.type_constraints;'Smoothing '];

%### add bottom track constraint

if exist('bvel','var')

 if sum(isfinite(bvel))>0 

  btweight=ps.velerr./bvels;

  [A2,A1,d,ubot,iubot]=lainbott(dbot.*wm,bvel,btweight,A2,A1,d,ps.botfac);

  if length(ubot)>0

   dr.zbot=z(iubot);

   dr.ubot=real(ubot(:,1));

   dr.vbot=imagnan(ubot(:,1));

   dr.uerrbot=ubot(:,2);

   % make length of array unique

   while length(dr.zbot)==nt | length(dr.zbot)==nz 

    dr.zbot=[dr.zbot;nan];

    dr.ubot=[dr.ubot;nan];

    dr.vbot=[dr.vbot;nan];

    dr.uerrbot=[dr.uerrbot;nan];

   end 

  else

   psbot=0;

  end

  if ps.botfac>0

   disp([' weight for bottom track is (ps.botfac) ',num2str(ps.botfac)])

   psbot=1;

  else

   disp(' not enought bottom track data')

   psbot=0; 

  end

 else

  disp(' no bottom track data')

  psbot=0; 

 end

else

 psbot=0;

end

de.ocean_constraints=[de.ocean_constraints;sum(abs(A2))-sum(de.ocean_constraints)];

de.ctd_constraints=[de.ctd_constraints;sum(abs(A1))-sum(de.ctd_constraints)];

de.type_constraints=[de.type_constraints;'Bottomtrk '];

%### add ship adcp constraint

if existf(di,'svel')==1

 if sum(isfinite(di.svel(:,1)))>2 

  [p,A2,A1,d,ds]=lainsadcp(p,di.svel,A2,A1,d,ps.dz,ps.sadcpfac,ps.velerr);

  if length(ds.z_sadcp)>1

   dr.z_sadcp=ds.z_sadcp;

   dr.u_sadcp=ds.u_sadcp;

   dr.v_sadcp=ds.v_sadcp;

   dr.uerr_sadcp=ds.uerr_sadcp;

   % make length of array unique

   if existf(dr,'zbot'), lzbot=length(dr.zbot); else, lzbot=0; end

   while length(dr.z_sadcp)==nt | length(dr.z_sadcp)==nz | ...

         length(dr.z_sadcp)==lzbot

    dr.z_sadcp=[dr.z_sadcp;nan];

    dr.u_sadcp=[dr.u_sadcp;nan];

    dr.v_sadcp=[dr.v_sadcp;nan];

    dr.uerr_sadcp=[dr.uerr_sadcp;nan];

   end

  end

  if ps.sadcpfac>0

   disp([' weight for SADCP vel is (ps.sadcpfac) ',num2str(ps.sadcpfac)])

  else

   disp(' no sadcp used')

   ps.sadcpfac=0;

  end

 else

  disp(' not enough SADCP data')

  ps.sadcpfac=0;

 end

else

  disp(' no SADCP data')

  ps.sadcpfac=0;

end

de.ocean_constraints=[de.ocean_constraints;sum(abs(A2))-sum(de.ocean_constraints)];

de.ctd_constraints=[de.ctd_constraints;sum(abs(A1))-sum(de.ctd_constraints)];

de.type_constraints=[de.type_constraints;'Ship ADCP '];

%### add barotropic constraint 

% check if position data exist 

uship_a=p.uship+sqrt(-1)*p.vship;

if (abs(uship_a)==0 & p.lat==0 & p.lon==0) | ~isfinite(p.lon+p.lat)

  disp(' no position data ');

  ps.barofac=0;

  ps.dragfac=0;

end

if ps.barofac>0

 fac=ps.velerr/ps.barvelerr;

 % need to increase if parts of profile are missing

 ii=find(di.dt>3*mean(di.dt));

 if length(ii)>1

  facgap=sum(di.dt(ii))/sum(di.dt);

  disp([' lainbaro: ',int2str(facgap*100),'% of profile have no useful data '])

  fac=fac*(1-tanh(facgap/0.15));

 end

 if ~isfinite(fac), fac=1; end

 [A2,A1,d]=lainbaro(A2,A1,d,uship_a,dtiv,ps.barofac*fac);

end

de.ocean_constraints=[de.ocean_constraints;sum(abs(A2))-sum(de.ocean_constraints)];

de.ctd_constraints=[de.ctd_constraints;sum(abs(A1))-sum(de.ctd_constraints)];

de.type_constraints=[de.type_constraints;'GPS naviga'];

%### small deep ocean velocity

if sum(smallfac(:,2))>0

 [A2,A1,d]=lainsmal(A2,A1,d,smallfac);

end

de.ocean_constraints=[de.ocean_constraints;sum(abs(A2))-sum(de.ocean_constraints)];

de.ctd_constraints=[de.ctd_constraints;sum(abs(A1))-sum(de.ctd_constraints)];

de.type_constraints=[de.type_constraints;'Small flow'];

%### check if problem is well constrained

if (psbot==0 & ps.barofac==0 & ps.sadcpfac==0), 

    disp(' no bottom no barotropic no SADCP constrain => will set mean U,V to zero')

    dr.onlyshear=1;

    [A2,A1,d]=lainocean(A2,A1,d);

elseif existf(dr,'onlyshear') 

    dr=rmfield(dr,'onlyshear');

end

de.ocean_constraints=[de.ocean_constraints;sum(abs(A2))-sum(de.ocean_constraints)];

de.ctd_constraints=[de.ctd_constraints;sum(abs(A1))-sum(de.ctd_constraints)];

de.type_constraints=[de.type_constraints;'Zero mean '];

%### add CTD drag constraint

if ps.dragfac>0

 disp([' weight for drag is (ps.dragfac) ',num2str(ps.dragfac)])

 [A2,A1,d]=laindrag(A2,A1,d,ctdvel,ps.dragfac);

end

de.ocean_constraints=[de.ocean_constraints;sum(abs(A2))-sum(de.ocean_constraints)];

de.ctd_constraints=[de.ctd_constraints;sum(abs(A1))-sum(de.ctd_constraints)];

de.type_constraints=[de.type_constraints;'Drag Model'];

[ld,a1l]=size(A1);

[ld,a2l]=size(A2);

disp([' ready for inversion  length of  d: ',num3str(ld,6,0)])

disp(['           (CTD vel)  length of A1: ',num3str(a1l,6,0)])

disp(['         (ocean vel)  length of A2: ',num3str(a2l,6,0)])