% Compute plasma current distribution Iy(R,Z)
% Inputs = {Fx(R,Z), ipsign, ag}
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function [Iy, varargout] = compute_Iy(Fx, ipsign, L, ag)
%% Axis and LCFS flux
[rA,zA,FA,dr2FA,dz2FA,drzFA,rX,zX,FX,dr2FX,dz2FX,drzFX,rl,zl,...
          rB,zB,FB,lB,lX,Opy] = meqpdom(Fx,ipsign,L.P.isaddl,L);

if isnan(FB)
    % Safe exit for optimization routine
    Iy = NaN(L.nzy, L.nry);
    varargout{1} = NaN;
    return
end
if rB == max(rl)
   % reject solutions limited touching the outer limiter
   Iy = NaN(L.nzy, L.nry);
   varargout{1} = NaN;
    return
end

%% Plasma current from basis function coefficients
[Tyg,~,~] = L.P.bfct(1,L.P.bfp,Fx,FA,FB,Opy,L.ry,L.iry);
Iy = reshape(Tyg*ag,L.nzy,L.nry);

%% Remove current outside boundary
Iy(~Opy) = 0;

%% Add outputs
if nargout >1
   % I woul like to use meqlarg for this purpose but path problem due to fact that it is inside a `private` folder
   out = struct('rA', rA, 'zA', zA, 'FA', FA, 'FB', FB, 'rX', rX, 'zX', zX, 'FX', FX,...
                'rB', rB, 'zB', zB, 'zl', zl, 'rl', rl, 'lB', lB, 'lX', lX, 'Opy', Opy,...
                 'Ip', sum(sum(Iy)), 'Tyg', Tyg, 'dr2FA', dr2FA, 'dz2FA', dz2FA, 'drzFA', drzFA);
   varargout{1} = out;
end
return