function [dFydag,dFprzdag,dFdIp,dFprzdIe,dFprzdIp,RBr,RBz,Gamma0] = rzpFlinfast(L,LY)
%% [dFdIy,dFdag,dFdIp,dFdIe,dFdCo,RBr,RBz,Gamma0] = rzpFlinfast(L,LY)
% Fast version of the rzip linear force balance calculation. The reduced
% force balance is F=[Fy;Fp;Fr;Fz] for the linearisation wrt ag (but dFydag
% is stored separately from the rest of the linearisation to ease the
% matrix inversion in rzpfastA for Xee calculation), F=[Fp;Fr;Fz] for Ie,
% while for the linearisation wrt Co only Fp is used. For a complete 
% discussion about the linearisation of the rzip force balance, refer to 
% rzpFlin. 
%
% This function calculates only the necessary elements for the
% calculation of the matrix A, and is not supposed to be used in a full
% state space calculation. 
% 
% It is called by rzpfastA and mexified as
% well.
%
% Input
%   L: L structure from rzp
%   LY: LY structure from LY, single time slice
%
% Output
%   dFydag: contains only dIydag = [dIydIp,dIydRc,dIydZc]
%   dFdag: linearisation of Fp = dIp - dIp', Fr and Fz wrt ag = [Ip,Rc,Zc]
%   dFdIp: placeholder for linearisation wrt Ip in case of a cde (never used)
%   dFdIe: linearisation of Fp, Fr, Fz wrt Ie
%   dFdCo: contains only the linearisation of Fp wrt Cp = [Ip,bp,li]
%   RBr,RBz,Gamma0: stores 2*pi*R*Br adn 2*pi*R*Bz and Gamma0. Shorten 
%                   computation time in rzpfastA 
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%% Store linearisation points
%#codegen
Ia = LY.Ia;
Iu = LY.Iu;
Iy = LY.Iy;
Ip = LY.Ip;
rIp = LY.rIp;
RBrye = -L.RBrye;
RBzye =  L.RBzye;
dr = L.drx;
dz = L.dzx;
ny = L.ny;
np = L.np;
ne = L.ne;

%% Compute analytical quantities of interest
%calculate the spatial derivatives of Iy with respect to R and Z
%with second order finite difference 
[dIydr, dIydz] = meqIyfd(Iy,dr,dz);

%Store Ie
Ie = [Ia;Iu];

RBz = RBzye*Ie;
RBr = RBrye*Ie;

%Gamma0 definition to obtain the equilibrium Fhoop = - Fr,ext
Gamma0 = -4*pi*( Iy(:)'*RBz(:) )/(mu0*Ip^2); 
      
%% Compute fast version linearization w.r.t Iy
dIydIp =  Iy(:)./Ip;

dIydRc = -dIydr(:);
dIydZc = -dIydz(:);        

dFydag = -[dIydIp dIydRc dIydZc];

%% Compute fast version linearisations w.r.t ag
dFrdIp = 2*pi*( Iy(:)' *RBz )/Ip + mu0*Ip*Gamma0;
dFrdRc = 2*pi*( dIydRc'*RBz )    + mu0*(Ip.^3)/(2*rIp);
dFrdZc = 2*pi*( dIydZc'*RBz );
dFzdRc = 2*pi*( dIydRc'*RBr );
dFzdZc = 2*pi*( dIydZc'*RBr );

dFprzdag = [     1     0       0;...
            dFrdIp dFrdRc dFrdZc;...
                0  dFzdRc dFzdZc]; 

%% Compute fast version linearisation w.r.t Ie
dFprzdIe = [zeros(1,ne);
            2*pi*(Iy(:)'*RBzye);...
            2*pi*(Iy(:)'*RBrye);];
      
%% Compute fast version linearisation w.r.t Co
dFprzdIp = [-1; 0; 0]; 

dFdIp = zeros(ny,np); %placeholder, still not implemented for the fast calculation

end