function [L] = rzpl(L,LX)
% RZP linearization around initial equilibrium guess LX.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

assert(nargin==2,'not enough input arguments');

L = fgel(L,LX);

lin = L.lin;

%% Store linearization points
lin.ag = [LX.Ip,LX.rIp./LX.Ip,LX.zIp./LX.Ip]; %expansion parameters for Iy in rzip are [Ip,Rc,Zc]

Bz = [L.G.Bzxa(L.lxy(:),:) L.G.Bzxu(L.lxy(:),:)]*lin.Ie;
Gamma0 = -4*pi*L.rry(:)'.*LX.Iy(:)'*Bz(:)/(mu0*lin.Ip^2);
%definition from A. Coutlis, Nucl. Fusion, vol. 39, no. 5, pp. 663–683, 1999.
%Linked to the same paper definition of Gamma0 in rzpFlin
lin.Lp = mu0*(lin.rIp/lin.Ip)*(Gamma0-0.5-lin.bp);

%% constraining quantities in LX

if all(L.icde) 
  switch L.P.cde
    case 'OhmTor_rigid'
      dLpdr = mu0*(Gamma0-0.5-LX.bp + 1); dLpdz = 0; %as in FP thesis
      dcdeIpdot = lin.Lp; %linear cde wrt Ip
      dcdeRcdot = lin.Ip*dLpdr + 2*pi*sum(L.rry(:).*LX.Iy(:).*Bz(:))/(lin.Ip); %linear cde wrt Rc
      dcdeZcdot = lin.Ip*dLpdz; %linear cde wrt Zc
      dcdedag = [dcdeIpdot dcdeRcdot dcdeZcdot];
      lin.Xpe = dcdedag * lin.dagdIe;
      lin.Xpo = dcdedag * lin.dagdCo;
    otherwise
      % placehold for implementation of different CDE linearization
  end
end

L.lin = lin;


