function [dIydIe,dIydCo,dFxdIe,dFxdCo] = fgslin(L,LY)
% function [dIydIe,dIydCo] = fgslin(L,LY)
% Returns linearization of Iy w.r.t. FGS inputs Ia, Iu, c(onstraints)  
% Optional output arguments [dIydIe,dIydCo,dFxdIe,dFxdCo] = fgslin(L,LY)
% also computes derivatives of Fx
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%% get derivatives of operator
[dFdIy,dFdag,dFdIp,dFdIe,dFdCo] = fgsFlin(L,LY);
[dIydIe,dIydCo,dagdIe,dagdCo] = fgselim(L,dFdIy,dFdag,dFdIp,dFdIe,dFdCo);

%%
if nargout>2
  %% get dFx/dIa by solving GS equation for perturbed dIa, dIy/dIa*dIa
  dFxdIe = zeros(L.nx,L.ne);
  for ie = 1:L.ne
    dIy = reshape(dIydIe(:,ie),L.nzy,L.nry); II = eye(L.ne); dIei=II(:,ie);
    dFxdIe(:,ie) = reshape(meqFx(L,dIy,dIei),[],1);
  end
  
  % get dFx/dCo
  dFxdCo = zeros(L.nx,L.nC);
  for io = 1:L.nC
    dIy = reshape(dIydCo(:,io),L.nzy,L.nry);
    dFxdCo(:,io) = reshape(meqFx(L,dIy,zeros(L.ne,1)),[],1);
  end
end
end
