function A = rzpfastA(Ia,Iu,Iy,Ip,rIp,bp,Mey,Mee,Re,Rp,RBrye,RBzye,drx,dzx,ny,np,ne,icde)
%% A = rzpfastA(Ia,Iu,Iy,Ip,rIp,bp,Mey,Mee,Re,Rp,RBrye,RBzye,rry,zzy,ny,np,ne,icde)
% Function to be mexified. It comprises the fast version for the system
% matrix (A) calculation. The matrix A calculation is the equivalent to the
% one performed in fgess, however it only uses the necessary steps,
% speeding up the computation.
%
% An example of the use of the function and how to retrieve the matrix A
% from fgess follow:
%
% [L,LX] = rzp(tok,shot,t);
%
% %% Using rzpfastA
% A = rzpfastA(LX.Ia,LX.Iu,LX.Iy,LX.Ip,LX.rIp,LX.bp,...
%         L.Mey,L.Mee,L.Re,LX.Rp,L.RBrye,L.RBzye,L.drx,L.dzx,L.ny,L.np,L.ne,L.icde);
%
% %% Using fgess
% sys = fgess(L);
% A = sys.A;
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%#codegen  
L  = struct('drx',drx,'dzx',dzx,'ny',ny,'np',np,'ne',ne,'RBrye',RBrye,'RBzye',RBzye);
LY = struct('Ia',Ia,'Iu',Iu,'Iy',Iy,'Ip',Ip,'rIp',rIp);

%% Force balance linearisation
[dFydag,dFprzdag,~,dFprzdIe,dFprzdIp,~,RBz,Gamma0] = rzpFlinfast(L,LY);

%% Compute derivatives of ag wrt Ie (short version of fgselim)
dagdIe = dFprzdag\dFprzdIe;

%% Green's function matrix modification terms via plasma influence in circuit eq
% 0  = (Mee+Xee)*Idot     + Xeo * codot + RIe - Ve 
dcedag = Mey*dFydag; % derivative od the circuit eq wrt ag
Xee = dcedag*dagdIe;

%% Define state related matrices S,K
S_ = Mee + Xee;
K_ = -diag(Re);

%% CDE OhmtTor_rigid related quantities
if icde
  dagdIp = (dFprzdag\dFprzdIp);    %Derivatives of ag wrt Ip
  
  %Green's function matrix modification due to total plasma current in circuit eq
  Xep = dcedag*dagdIp;

  %cde linearisation wrt to ag
  dcdeIpdot = mu0*rIp*(Gamma0-0.5-bp)/Ip; %Lp0
  dcdeRcdot = mu0* Ip*(Gamma0-0.5-bp+1) + 2*pi*(Iy(:)'*RBz(:))/(Ip); %dLp0dRc + 2*pi*(Iy(:)'*RBz(:))/(Ip)
  dcdedag = [dcdeIpdot dcdeRcdot 0]; %dLp0dZc = 0
  
  %Green's function matrix modification due to plasma and total plasma current in cde
  Xpe = -dcdedag * dagdIe;
  Xpp = -dcdedag * dagdIp;

  %Add Ip state to state space 
  S = [S_, Xep; -dcedag(:,1)' + Xpe, Xpp]; % -dcedag(:,1)' = Mpe
  K = [K_, zeros(ne,1); zeros(1,ne), -Rp]; 
else
  S = S_;
  K = K_;
end

%% Generate A matrix in continuous time representation
A = S\K;

end