function [dIydIe,dIydCo,dagdIe,dagdCo] = fgselim(L,dFdIy,dFdag,dFdIp,dFdIe,dFdCo)
%% Compute derivatives of FGS solutions (Iy,ag) with respect to inputs (Ie,Co), in static case (no CDE), from derivatives of the operator fgsF 
%
% When used within FGE, rows of the operator linked to the circuit equations or the CDE are ignored.
% If the CDE is used, Ip needs to be treated as a constraint input rather than part of the state.
%
% Formulation:
% Given the FGS operator F = F(Iy,ag;Ie,Co)
% Solutions of the problem, for given inputs (Ie,Co), are formally Iy = Iy(Ie,co), ag = ag(Ie,co), such that F(Iy,ag;Ie,Co) = 0
%
% Linear approximation of the solution for perturbed inputs dIe,dCo around linearization point Ie0,Co0 reads
% dIy(dIe,dCo) = Iy(Ie0 + dIe,Co0 + dCo) - Iy(Ie0,Co0) = dIydIe * dIe + dIydCo * dCo
% dag(dIe,dCo) = ag(Ie0 + dIe,Co0 + dCo) - ag(Ie0,Co0) = dagdIe * dIe + dagdCo * dCo
%
% Identify dIydIe, dIydCo, dagdIe, dagdCo from derivatives of the operator F:
% a) compute the linear variation of the operator F
% [dFdIy,dFdag,dFdIe,dFdCo]*[dIy;dag;dIe;dCo] = 0, 
% where F(Iy0,ag0,Ie0,co0) = 0, and F(Iy+dIy,ag+dag,Ie+dIe,Co+dICo) = 0
%
% b) Solve for [dIy, dag],  bringing to the RHS dIe,dCo
% [dFdIy,dFdag]*[dIy;dag] = - [dFdIe, dFdCo]*[dIe; dCo]
%
% c) Calling M = -[dFdIy,dFdag]\[dFdIe, dFdCo],
% one can identify for example dIydIe = M(1:ny,1:ne), dagdIe = M(ny+1:end,1:ne)
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

iy = L.ind.iy; ig = L.ind.ig; 
je = 1:L.G.na + L.G.nu; jco=je(end)+(1:L.nC);

% In case of CDE, Ip perturbation occurs through Ip state
if L.icde
  iagconIp = contains(L.agconc(:,3),'Ip');
  assert(~any(dFdCo(:,iagconIp)),'expected dFdCo to be only zero here..');
  dFdCo(:,iagconIp) = dFdIp; % replace with dFdIp value
end


M = -[dFdIy([iy,ig],:), dFdag([iy,ig],:)]\[dFdIe([iy,ig],:), dFdCo([iy,ig],:)];

dIydIe = M(1:L.ny,je     );
dIydCo = M(1:L.ny,jco    );
dagdIe = M(L.ny+1:end,je ); 
dagdCo = M(L.ny+1:end,jco); 

end


