function [L] = fgel(L,LX)
% FGE linearization around initial equilibrium guess LX.
% Takes linearizations around the GS equilibrium using fgsFlin
% and computes linearizations for the coil current equation.
% and fills L.lin structure with linearization data to be used by fgess() 
% or fget.m
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

assert(nargin==2,'not enough input arguments');
assert(numel(LX.t)==1,'fgel works for one time slice only, use meqxk() to slice')

%if no plasma, no correction to inductance terms
if all(LX.Ip == 0)
  L.lin.Xee = zeros(L.ne);
  return
end

%% Compute linearization
if L.P.debug,  fprintf('Linearizing around equilibrium... '); end

%% fge/fgs/rzp Force operator Jacobians - dFdIp will be empty depending on use of CDE

[dFdIy,dFdag,dFdIp,dFdIe,dFdCo] = L.codeF(L,meqxk(LX,1));

lin.dFdIy = dFdIy; lin.dFdIe = dFdIe; lin.dFdIp = dFdIp;
lin.dFdag = dFdag; lin.dFdCo = dFdCo; 

%% Store linearization points
lin.ag = LX.ag;
lin.Ie = [LX.Ia;LX.Iu];
lin.Ia = LX.Ia;
lin.Iu = LX.Iu;
lin.Bm = LX.Bm;
lin.Ff = LX.Ff;
lin.rA = LX.rA;
lin.zA = LX.zA;
lin.zIp = LX.zIp;
lin.rIp = LX.rIp; 
lin.zIp = sum(L.zzy(:).*LX.Iy(:));
lin.rIp = sum(L.rry(:).*LX.Iy(:));
lin.Ip = LX.Ip;
if L.nn
  lin.Fn = LX.Fn; 
  lin.Brn = LX.Brn;
  lin.Bzn = LX.Bzn;
end
lin.Ini = LX.Ini;
lin.SQ  = LX.SQ;
lin.Fx  = LX.Fx; 
lin.FA  = LX.FA;
lin.FB  = LX.FB;
lin.Ft  = LX.Ft;
lin.bp  = LX.bp;
lin.qA  = LX.qA;
lin.li  = LX.li;

lin.Mpe = (L.Mey*LX.Iy(:))' ./ LX.Ip;
lin.Rp = LX.Rp;
lin.Lp = LX.Lp;


%% constraining quantities in LX
lin.Co = LX2Co(L,LX); %constraint quantities [Ip,bp,li/qa] (only li for rzip for now)

[dIydIe,dIydCo,dagdIe,dagdCo] = fgselim(L,dFdIy,dFdag,dFdIp,dFdIe,dFdCo);

% Green's function matrix modification terms via plasma influence in circuit eq
% 0  = (Mee+Xee)*Idot     + Xeo * codot + RIe - Ve
lin.Xee = L.Mey * dIydIe; lin.Xeo = L.Mey * dIydCo;

if all(L.icde) 
  lin.Xpe = zeros(L.np,L.ne);
  lin.Xpo = zeros(L.np,L.nC);
  switch L.P.cde
    case 'OhmTor_rigid'
      lin.Xpo(contains(L.agconc(:,3),'Ip')) = lin.Lp;
    otherwise
      % placehold for implementation of different CDE linearization
  end
end

% Bm = ([Bma,Bmu]+Xme)*(Ie-Ie0) + Xmo*(Co-Co0) + Bm0  % Fields
% Ff = ([Mfa,Mfu]+Xfe)*(Ie-Ie0) + Xfo*(Co-Co0) + Ff0  % Fluxes
lin.dBmdIe = [L.G.Bma,L.G.Bmu] + L.G.Bmx(:,L.lxy) * dIydIe; % Bm
lin.dFfdIe = [L.G.Mfa,L.G.Mfu] + L.G.Mfx(:,L.lxy) * dIydIe; % Ff
lin.dBmdCo =                     L.G.Bmx(:,L.lxy) * dIydCo;
lin.dFfdCo =                     L.G.Mfx(:,L.lxy) * dIydCo;

lin.dIadIe = [eye(L.G.na),zeros(L.G.na,L.G.nu)]; % Ia  from the state
lin.dIudIe = [zeros(L.G.nu,L.G.na),eye(L.G.nu)]; % Iu  from the state

lin.dagdCo = dagdCo;
lin.dagdIe = dagdIe;
lin.dIydCo = dIydCo;

[dFxdIy, dFAdIy, dFBdIy, drAdIy, dzAdIy, dFndIy, dBrndIy, dBzndIy, dFtdIy, dIpdIy, dbpdIy, dbpDdIy, dlidIy, dqAdIy,...
 dFxdIe, dFAdIe, dFBdIe ,drAdIe, dzAdIe, dFndIe, dBrndIe, dBzndIe, dFtdIe, dIpdIe, dbpdIe, dbpDdIe, dlidIe, dqAdIe,...
                                                                   dFtdag, dIpdag, dbpdag, dbpDdag, dlidag, dqAdag] = dFdI(L,LX);

dFABdIy = (dFAdIy + dFBdIy(1:LX.nA,:))/2;
dFABdIe = (dFAdIe + dFBdIe(1:LX.nA,:))/2;

lin.dFABdIe = dFABdIy*dIydIe + dFABdIe;
lin.dFABdCo = dFABdIy*dIydCo;
lin.dFAdIe  = dFAdIy*dIydIe + dFAdIe;
lin.dFAdCo  = dFAdIy*dIydCo;
lin.dFBdIe  = dFBdIy*dIydIe + dFBdIe;
lin.dFBdCo  = dFBdIy*dIydCo;

% Post-processing quantity linearizations
lin.dIydIe   = dIydIe;
lin.drAdIe   = drAdIe           + drAdIy *(dIydIe)          ;
lin.dzAdIe   = dzAdIe           + dzAdIy *(dIydIe)          ;
lin.dFndIe   = dFndIe           + dFndIy *(dIydIe)          ;
lin.dBrndIe  = dBrndIe          + dBrndIy*(dIydIe)          ;
lin.dBzndIe  = dBzndIe          + dBzndIy*(dIydIe)          ;
lin.dFxdIe   = dFxdIe           + dFxdIy *(dIydIe)          ;
lin.dFtdIe   = dFtdag *(dagdIe) + dFtdIy *(dIydIe) + dFtdIe ; 
lin.dlidIe   = dlidag *(dagdIe) + dlidIy *(dIydIe) + dlidIe ; 
lin.dbpdIe   = dbpdag *(dagdIe) + dbpdIy *(dIydIe) + dbpdIe ; 
lin.dbpDdIe  = dbpDdag*(dagdIe) + dbpDdIy*(dIydIe) + dbpDdIe; 
lin.dqAdIe   = dqAdag *(dagdIe) + dqAdIy *(dIydIe) + dqAdIe ; 
lin.dFxdCo   = dFxdIy *(dIydCo)                             ;
lin.dFtdCo   = dFtdIy *(dIydCo) + dFtdag *(dagdCo)          ;
lin.dbpdCo   = dbpdIy *(dIydCo) + dbpdag *(dagdCo)          ;
lin.dbpDdCo  = dbpDdIy*(dIydCo) + dbpDdag*(dagdCo)          ;
lin.dlidCo   = dlidIy *(dIydCo) + dlidag *(dagdCo)          ;
lin.dqAdCo   = dqAdIy *(dIydCo) + dqAdag *(dagdCo)          ;


% rIp, zIp, Ip variations due to Ie, Co: directly from dIydIe / dIydCo
lin.drIpdIe = sum(L.rry(:).*dIydIe,1);
lin.dzIpdIe = sum(L.zzy(:).*dIydIe,1);
lin.dIpdIe  = sum(dIydIe,1);
lin.drIpdCo = sum(L.rry(:).*dIydCo,1);
lin.dzIpdCo = sum(L.zzy(:).*dIydCo,1);
lin.dIpdCo  = sum(dIydCo,1);

% dbpdIe,dqAdIe must be 0 when bp,and qA are constraints.
% dbpdCo,dqAdCo must be 1 for the corresponding index when bp,and qA are constraints.
% From the linearization dbpdIe,dqAdIe~0, dbpdCo,dqAdCo~1
% but the small error may accumulate for long simulations because its value 
% is being fed back to compute the time derivative in fgetkl(). 
% Hence, forced to be the analytic values.

for ii = 1:L.nC
  if strcmp(L.agconc{ii,3},'Ip')
    continue
  end
  dCoiidIe = sprintf('d%sdIe',L.agconc{ii,3});
  dCoiidCo = sprintf('d%sdCo',L.agconc{ii,3});
  lin.(dCoiidIe) = zeros(1,L.ne); 
  lin.(dCoiidCo) = zeros(1,L.nC);
  lin.(dCoiidCo)(1,ii) = 1;
end

lin.dFndCo  = dFndIy*(dIydCo);
lin.dBrndCo = dBrndIy*(dIydCo);
lin.dBzndCo = dBzndIy*(dIydCo);

lin.drAdCo  = drAdIy*(dIydCo);
lin.dzAdCo  = dzAdIy*(dIydCo);

%% (TCV) doublet observers
if L.P.idoublet
  for ii=1:L.nD
    O = (LX.Opy(:)==ii);
    lin.drIpDdIe(ii,:) = sum(L.rry(:).*O.*dIydIe);
    lin.dzIpDdIe(ii,:) = sum(L.zzy(:).*O.*dIydIe);
    lin. dIpDdIe(ii,:) = sum(O.*dIydIe);
    lin.drIpDdCo(ii,:) = sum(L.rry(:).*O.*dIydCo);
    lin.dzIpDdCo(ii,:) = sum(L.zzy(:).*O.*dIydCo);
    lin. dIpDdCo(ii,:) = sum(O.*dIydCo);
  end
  lin.IpD = LX.IpD;
  lin.rIpD = LX.rIpD;
  lin.zIpD = LX.zIpD;
  lin.bpD = LX.bpD;
end

%%

L.lin = lin;

if L.P.debug,  fprintf('done\n'); end


