function [L,LX,LY] = lih(tok,shot,t,varargin)
%LIH Finite-element reconstruction of plasma current density
% Essentially based on first step of LIU
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% specific lih parameters different than liu - put them before any possible
% user-defined overrides.
PP = {'rop',NaN,'rip',NaN,'zup',NaN,'zlp',NaN,'nelem',10};

if nargout==1
  L = liu(tok,shot,t,PP{:},varargin{:});
elseif nargout>=2
  [L,LX] = liu(tok,shot,t,PP{:},varargin{:});
end
L = lihc(L); % custom lihc consolidation

if nargout==3
  % LIH reconstruction
  if L.P.slx
    LY = lihtsim(L,LX, PP{:},varargin{:}); % Simulink version
  else
    LY = liht(L,LX); % Matlab version
  end
end

end

function L=lihc(L)
% Additional consolidation for lihc (parts not already in liu)
L.code = 'lih';

% observers
mask_t = L.zzy(:) > 0; % mask for top    half of grid
mask_b = L.zzy(:) < 0; % mast for bottom half of grid

if any(L.zzy(:))
  % split contribution  - makes mask double
  mask_t = double(mask_t); mask_b=double(mask_b);
  mask_t(L.zzy(:)==0) = 1/2;
  mask_b(L.zzy(:)==0) = 1/2;
end

L.Iph    = ones(1,L.nh);
L.Ipzh   = sum(L.zzy(:)        .*L.Tyh,1);
L.Iprh   = sum(L.rry(:)        .*L.Tyh,1);
L.Iph_t  = sum(          mask_t.*L.Tyh,1);
L.Iph_b  = sum(          mask_b.*L.Tyh,1);
L.Ipzh_t = sum(L.zzy(:).*mask_t.*L.Tyh,1);
L.Ipzh_b = sum(L.zzy(:).*mask_b.*L.Tyh,1);
L.Iprh_t = sum(L.rry(:).*mask_t.*L.Tyh,1);
L.Iprh_b = sum(L.rry(:).*mask_b.*L.Tyh,1);
L.mask_t = mask_t;
L.mask_b = mask_b;

end
