function LY = liht(L,LX,varargin)
% LIUQE-like finite element plasma current distribution observer
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% optional parameters
for k = 1:2:length(varargin), L.P.(varargin{k}) = varargin{k+1}; end

nt = numel(LX.t);

for kt = 1:nt
  t  = LX.t(kt);
  Ip = LX.Ip(kt); rBt = LX.rBt(kt);
  shot = LX.shot(min(kt,end));
  
  Xr = [LX.Ff(:,kt); LX.Bm(:,kt)]; % magnetic measurements
  Xe = [LX.Ia(:,kt); LX.Iu(:,kt)]; % external currents
  Xi = [LX.Ft(kt)*rBt ; Ip];       % integral measurements
  Yr = L.Wr.*Xr;
  Ye = L.We.*Xe;
  Yi = L.Wi.*Xi;
  Yd = [Yr ; Ye ; Yi];
  
  if L.P.iterh
    Jh0 = repmat(Ip/L.nh,L.nh,1);
    [Ie,Jh,~] = ipmhmex(L.Ahd,L.Aed,L.Ahe,L.Aeh,L.Ahh,L.uAhh,Yd,Xe,Jh0,sign(Ip),L.P.iterh,max(L.P.Ipmin/L.nh,abs(Jh0(1)))*L.P.tolh);
  else
    Ie = L.Aed*Yd;
    Jh = L.Ahd*Yd;
  end
  Zd = L.Wdh*Jh+L.Wde*Ie; 
  chih = sqrt(sum((Zd-Yd).^2)/L.nd); % chi of current fit

  Ip      = L.Iph  * Jh ;
  zIp     = L.Ipzh * Jh ;
  rIp     = L.Iprh * Jh ;
  
  if L.P.idoublet
    IpD     = [L.Iph_t  * Jh ; L.Iph_b  * Jh];
    rIpD    = [L.Iprh_t * Jh ; L.Iprh_b * Jh];
    zIpD    = [L.Ipzh_t * Jh ; L.Ipzh_b * Jh];
  else
    IpD = Ip; zIpD=zIp; rIpD=rIp;
  end
  rA = rIpD./IpD;
  zA = zIpD./IpD;
  
  Ia = Ie(1:L.G.na);
  Iu = Ie(L.G.na+(1:L.G.nu));
  Iy = reshape(L.Tyh*Jh,L.nzy,L.nry);
  
  %% Poisson solver for a coarse flux estimate
  Fx = meqFx(L,Iy,Ie,0,1,Jh);
  if L.P.ivacuum, end

  %% mask of droplet viability
  nA = sum(abs(IpD) > L.P.Ipmin/2);
  
  %% Postprocessing, faking some plasma-related things
  TpDg = zeros(L.nD,L.ng); ITpDg=zeros(L.nD,L.ng); ag=zeros(L.ng,1); 
  [rB,zB,FB] = deal(zeros(L.nD,1));
  [rX,zX,FX] = deal([]);
  [FA,dr2FA,dz2FA,drzFA] = deal(zeros(nA,1));
  
  lB = true;
  [lX] = deal(false(L.nD,1));
  Opy = zeros(L.nzy,L.nry,'int8');
  
  LYt = struct('shot',shot,'t',t);
  [LYt] = meqpost(L,LYt,TpDg,ITpDg,ag,...
  Fx,FA,FB,rA,zA,dr2FA,dz2FA,drzFA,rB,zB,lB,lX,rX,zX,FX,...
  rBt,Ia,Iu,Iy,Opy);
  
  % override some values
  LYt.Ip = sum(Iy); 
  LYt.IpD = IpD;   
  LYt.rIpD = rIpD;
  LYt.zIpD = zIpD;
  LYt.nA = nA;      
  
  % lih specials
  LYt  = meqlarg(LYt,chih,Ip); 
  
  LY(kt) = LYt; %#ok<AGROW>
end

LY = meqlpack(LY);

% Time derivatives
if numel(LY.t)>1
  dt = diff(LY.t);
  LY.Iadot = [zeros(L.G.na,1)          ,diff(LY.Ia,[],2)./dt];
  LY.Iudot = [zeros(L.G.nu,1)          ,diff(LY.Iu,[],2)./dt];
  LY.Iydot = cat(3,zeros(L.nzy,L.nry,1),diff(LY.Iy,[],3)./permute(dt,[1 3 2]));
else
  LY.Iadot = zeros(L.G.na,1);
  LY.Iudot = zeros(L.G.nu,1);
  LY.Iydot = zeros(L.nzy,L.nry,1);
end
LY.Uf = L.G.Mfa*LY.Iadot + L.G.Mfu*LY.Iudot + L.G.Mfx(:,L.lxy(:))*reshape(LY.Iydot,[],nt);

end
