function Fx = meqFx(L,IyD,Ie,dz,rst,Jh)
%MEQFX Evaluate flux for given Iy, and optionally Ie
% FX = MEQFX(L,IYD[,IE,DZ,RST,JH]) returns the solution of the poisson
% equation for given current density parametrization. The L structure
% contains MEQ ancillary data, IYD contains the breakdown of the plasma
% current per plasma domain , IE is the current in the external conductors,
% assumed zero if not passed.
% The next 3 arguments are optional and are typically used only in LIU. DZ
% is the vector of shifts of IY with respect to FX for each domain, if RST
% is true then JH contains the current in each FE for the LIU initial
% guess.
%
% Note: outside of LIU or if L.ndz<=1 then IY can be used in place of IYD.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if nargin<3
  Iyonly = true; % compute only plasma contribution
else
  Iyonly = false; % include extenal conductors
end

if all(IyD(:)==0)
  Ieonly = true;
else
  Ieonly = false;
end

if nargin<6 || rst || ~any(dz)
  if nargin<6, rst = false;end
  assert(~(rst && any(dz)),'If rst is true, all dzs should be 0');
  dz = 0;
  dzMbe = zeros(1,L.ne);
  ndzeff = int8(1);
  IyD(:,:,1) = sum(IyD,3);
else
  ndzeff = find(any(any(IyD,1),2),1,'last');
  if isempty(ndzeff), ndzeff = 1; end % No plasma current
  switch L.P.gsxe
    case 1, dzMbe = L.dzMbe;
    case 2, dzMbe = L.dzMboe;
  end
end

assert(ndzeff <= numel(dz),'numel(dz)=%d, but ndzeff=%d',numel(dz),ndzeff);

% First domain
Iy = IyD(:,:,1);
if ~Ieonly
  if rst % from FE
    Fb = L.Mbh*Jh;
  else % from Iy, as sheet current equivalent to normal flux derivative
    Fb = meqfbp(Iy,L); % boundary current contribution from plasma - Lackner's trick!
  end
else
  Fb = zeros(L.nx-L.ny,1);
end

switch L.P.gsxe
  case 1 % Include external currents inside the grid in Iy - LIUQE mode
    if ~Iyonly
      Fb = Fb + L.Mbe*Ie - dzMbe*(Ie*dz(1)); % Add effect of external currents and dz
      Iy(L.kyw) = L.Tww*Ie(L.kew);        % Add external currents on computational grid
    end
    % Note: This does not take dz into account
    Fx = gszrmex(Fb,Iy,L.cx,L.cq,L.cr,L.cs,L.ci,L.co,L.idzx*dz(1));
  case 2 % Treat external currents inside the grid exactly - FBT mode
    if ~Iyonly
      Fb = Fb + L.Mboe*Ie(L.ioe) - dzMbe*(Ie*dz(1)); % Add effect of external currents and dz
    end
    Fx = gszrmex(Fb,Iy,L.cx,L.cq,L.cr,L.cs,L.ci,L.co,L.idzx*dz(1)); % Solve GS
    
    if ~isempty(L.Mxie) && ~Iyonly
      Fxie = L.Mxie*Ie(L.iie); % Flux from external currents in the grid
      Fx = Fx + reshape(Fxie,L.nzx,L.nrx); % add flux
    end
  otherwise
    error('meqFX:invalidmode','gsex=%d is invalid',L.P.gsex);
end

% Other domains (no external currents)
for iD = 2:ndzeff
  Iy = IyD(:,:,iD);
  Fb = meqfbp(Iy,L);
  Fx = Fx + gszrmex(Fb,Iy,L.cx,L.cq,L.cr,L.cs,L.ci,L.co,L.idzx*dz(iD)); % Solve GS
end
end