function [Fxy, FAy, FBy, rAy, zAy, Fny, Brny, Bzny, Fty, Ipy, bpy, bpDy, liy, qAy,...
          Fxe, FAe, FBe ,rAe, zAe, Fne, Brne, Bzne, Fte, Ipe, bpe, bpDe, lie, qAe,...
                                                    Ftg, Ipg, bpg, bpDg, lig, qAg,...
          Sy,Se,Sg] = dFdI(L,LY)
% Evaluate derivatives of post-processing outputs of F operator
% w.r.t. changes in Iy, Ia, Iu
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% Enable computation of signeo derivative. 
% Since it involves expensive meqpostq, done only when necessary
dsig = nargout > 34; 

assert(numel(LY.t)==1,'dFdI must be called for 1 timeslice only, use mexqk() to slice')

Ia = LY.Ia; Iu = LY.Iu;
Iy = LY.Iy;
rBt = LY.rBt; ag = LY.ag;
signeo = LY.signeo;

dI = 1e-3*max(abs(Iy(:)));

% init
Fxy = zeros(L.nx,    L.ny);
FAy = zeros(LY.nA, L.ny);
FBy = zeros(LY.nB, L.ny);
rAy = zeros(LY.nA, L.ny);
zAy = zeros(LY.nA, L.ny);
Fny = zeros(L.nn,    L.ny);
Brny = zeros(L.nn,    L.ny);
Bzny = zeros(L.nn,    L.ny);
Fty = zeros(1,L.ny);
Ipy = zeros(1,L.ny);
bpy = zeros(1,L.ny);
bpDy = zeros(L.nD,L.ny);
liy = zeros(1,L.ny);
qAy = zeros(LY.nA,L.ny);
Sy  = zeros(LY.nB, L.ny);
for ii=1:numel(Iy)
  dIy = zeros(L.nzy,L.nry);
  dIy(ii) = dI;
  
  if dsig
    [Fxp,FAp,FBp,rAp,zAp,Fnp,Brnp,Bznp,Ftp,Ipp,bpp,bpDp,lip,qAp,Sp] = evalFxy(L,Iy+dIy,Ia,Iu,rBt,ag,signeo);
    [Fxm,FAm,FBm,rAm,zAm,Fnm,Brnm,Bznm,Ftm,Ipm,bpm,bpDm,lim,qAm,Sm] = evalFxy(L,Iy-dIy,Ia,Iu,rBt,ag,signeo);
    Sy(:,ii)  = (Sp(:)-Sm(:))/(2*dIy(ii));
  else
    [Fxp,FAp,FBp,rAp,zAp,Fnp,Brnp,Bznp,Ftp,Ipp,bpp,bpDp,lip,qAp   ] = evalFxy(L,Iy+dIy,Ia,Iu,rBt,ag       );
    [Fxm,FAm,FBm,rAm,zAm,Fnm,Brnm,Bznm,Ftm,Ipm,bpm,bpDm,lim,qAm   ] = evalFxy(L,Iy-dIy,Ia,Iu,rBt,ag       );
  end
  Fxy(:,ii) = (Fxp(:)-Fxm(:))/(2*dIy(ii));
  FAy(:,ii) = (FAp-FAm)/(2*dIy(ii));
  FBy(:,ii) = (FBp-FBm)/(2*dIy(ii));
  rAy(:,ii) = (rAp-rAm)/(2*dIy(ii));
  zAy(:,ii) = (zAp-zAm)/(2*dIy(ii));
  Fny(:,ii) = (Fnp-Fnm)/(2*dIy(ii));
  Brny(:,ii)= (Brnp-Brnm)/(2*dIy(ii));
  Bzny(:,ii)= (Bznp-Bznm)/(2*dIy(ii));
  Fty(:,ii) = (Ftp-Ftm)/(2*dIy(ii));
  Ipy(:,ii) = (Ipp - Ipm)/(2*dIy(ii));
  bpy(:,ii) = (bpp - bpm)/(2*dIy(ii));
  bpDy(:,ii) = (bpDp - bpDm)/(2*dIy(ii));
  liy(:,ii) = (lip - lim)/(2*dIy(ii));
  qAy(:,ii) = (qAp - qAm)/(2*dIy(ii));
end

% w.r.t. changes to Ia
Fxa = zeros(L.nx, L.G.na);
FAa = zeros(LY.nA, L.G.na);
FBa = zeros(LY.nB, L.G.na);
rAa = zeros(LY.nA, L.G.na);
zAa = zeros(LY.nA, L.G.na);
Fna = zeros(L.nn, L.G.na);
Brna = zeros(L.nn, L.G.na);
Bzna = zeros(L.nn, L.G.na);
Fta = zeros(1, L.G.na);
Ipa = zeros(1,L.ng);
bpa = zeros(1,L.ng);
bpDa = zeros(L.nD,L.ng);
lia = zeros(1,L.ng);
qAa = zeros(LY.nA,L.ng);
Sa =  zeros(LY.nB, L.G.na);
for ii=1:L.G.na
  dIa = zeros(L.G.na,1);
  dIa(ii) = 1e-3;
  if dsig
    [Fxp,FAp,FBp,rAp,zAp,Fnp,Brnp,Bznp,Ftp,Ipp,bpp,bpDp,lip,qAp,Sp] = evalFxy(L,Iy,Ia+dIa,Iu,rBt,ag,signeo);
    [Fxm,FAm,FBm,rAm,zAm,Fnm,Brnm,Bznm,Ftm,Ipm,bpm,bpDm,lim,qAm,Sm] = evalFxy(L,Iy,Ia-dIa,Iu,rBt,ag,signeo);
    Sa(:,ii)  = (Sp(:)-Sm(:))/(2*dIa(ii));
  else
    [Fxp,FAp,FBp,rAp,zAp,Fnp,Brnp,Bznp,Ftp,Ipp,bpp,bpDp,lip,qAp  ] = evalFxy(L,Iy,Ia+dIa,Iu,rBt,ag       );
    [Fxm,FAm,FBm,rAm,zAm,Fnm,Brnm,Bznm,Ftm,Ipm,bpm,bpDm,lim,qAm  ] = evalFxy(L,Iy,Ia-dIa,Iu,rBt,ag       );
  end
  Fxa(:,ii) = (Fxp(:)-Fxm(:))/(2*dIa(ii));
  FAa(:,ii) = (FAp-FAm)/(2*dIa(ii));
  FBa(:,ii) = (FBp-FBm)/(2*dIa(ii));
  rAa(:,ii) = (rAp-rAm)/(2*dIa(ii));
  zAa(:,ii) = (zAp-zAm)/(2*dIa(ii));
  Fna(:,ii) = (Fnp-Fnm)/(2*dIa(ii));
  Brna(:,ii) = (Brnp-Brnm)/(2*dIa(ii));
  Bzna(:,ii) = (Bznp-Bznm)/(2*dIa(ii));
  Fta(:,ii) = (Ftp-Ftm)/(2*dIa(ii));
  Ipa(:,ii) = (Ipp - Ipm)/(2*dIa(ii));
  bpa(:,ii) = (bpp - bpm)/(2*dIa(ii));
  bpDa(:,ii) = (bpDp - bpDm)/(2*dIa(ii));
  lia(:,ii) = (lip - lim)/(2*dIa(ii));
  qAa(:,ii) = (qAp - qAm)/(2*dIa(ii));
end

% w.r.t. changes in Iu
Fxu = zeros(L.nx, L.G.nu);
FAu = zeros(LY.nA, L.G.nu);
FBu = zeros(LY.nB, L.G.nu);
rAu = zeros(LY.nA, L.G.nu);
zAu = zeros(LY.nA, L.G.nu);
Fnu = zeros(L.nn, L.G.nu);
Brnu = zeros(L.nn, L.G.nu);
Bznu = zeros(L.nn, L.G.nu);
Ftu = zeros(1, L.G.nu);
Ipu = zeros(1,L.G.nu);
bpu = zeros(1,L.G.nu);
bpDu = zeros(L.nD,L.G.nu);
liu = zeros(1,L.G.nu);
qAu = zeros(LY.nA,L.G.nu);
Su  = zeros(LY.nB, L.G.nu);
for ii=1:L.G.nu
  dIu = zeros(L.G.nu,1);
  dIu(ii) = 1;
  if dsig
    [Fxp,FAp,FBp,rAp,zAp,Fnp,Brnp,Bznp,Ftp,Ipp,bpp,bpDp,lip,qAp,Sp] = evalFxy(L,Iy,Ia,Iu+dIu,rBt,ag,signeo);
    [Fxm,FAm,FBm,rAm,zAm,Fnm,Brnm,Bznm,Ftm,Ipm,bpm,bpDm,lim,qAm,Sm] = evalFxy(L,Iy,Ia,Iu-dIu,rBt,ag,signeo);
    Su(:,ii)  = (Sp(:)-Sm(:))/(2*dIu(ii));
  else
    [Fxp,FAp,FBp,rAp,zAp,Fnp,Brnp,Bznp,Ftp,Ipp,bpp,bpDp,lip,qAp   ] = evalFxy(L,Iy,Ia,Iu+dIu,rBt,ag       );
    [Fxm,FAm,FBm,rAm,zAm,Fnm,Brnm,Bznm,Ftm,Ipm,bpm,bpDm,lim,qAm   ] = evalFxy(L,Iy,Ia,Iu-dIu,rBt,ag       );
  end
  Fxu(:,ii) = (Fxp(:)-Fxm(:))/(2*dIu(ii));
  FAu(:,ii) = (FAp-FAm)/(2*dIu(ii));
  FBu(:,ii) = (FBp-FBm)/(2*dIu(ii));
  rAu(:,ii) = (rAp-rAm)/(2*dIu(ii));
  zAu(:,ii) = (zAp-zAm)/(2*dIu(ii));
  Fnu(:,ii) = (Fnp-Fnm)/(2*dIu(ii));
  Brnu(:,ii) = (Brnp-Brnm)/(2*dIu(ii));
  Bznu(:,ii) = (Bznp-Bznm)/(2*dIu(ii));
  Ftu(:,ii) = (Ftp-Ftm)/(2*dIu(ii));
  Ipu(:,ii) = (Ipp - Ipm)/(2*dIu(ii));
  bpu(:,ii) = (bpp - bpm)/(2*dIu(ii));
  bpDu(:,ii) = (bpDp - bpDm)/(2*dIu(ii));
  liu(:,ii) = (lip - lim)/(2*dIu(ii));
  qAu(:,ii) = (qAp - qAm)/(2*dIu(ii));
end

% combine
Fxe = [Fxa,Fxu];
FAe = [FAa,FAu];
FBe = [FBa,FBu];
rAe = [rAa,rAu];
zAe = [zAa,zAu];
Fne = [Fna,Fnu];
Brne = [Brna,Brnu];
Bzne = [Bzna,Bznu];
Fte = [Fta,Ftu];
Ipe = [Ipa,Ipu];
bpe = [bpa,bpu];
bpDe = [bpDa,bpDu];
lie = [lia,liu];
qAe = [qAa,qAu];
Se  = [Sa,Su];



Ftg = zeros(1,L.ng);
Ipg = zeros(1,L.ng);
bpg = zeros(1,L.ng);
bpDg = zeros(L.nD,L.ng);
lig = zeros(1,L.ng);
qAg = zeros(LY.nA,L.ng);
Sg = zeros(LY.nB,L.ng);
dd = 1;
% w.r.t. changes in ag
for ii=1:L.ng
  dag = zeros(L.ng,1);
  dag(ii) = dd;
  if dsig
    [~,~,~,~,~,~,~,~,Ftp,Ipp,bpp,bpDp,lip,qAp,Sp] = evalFxy(L,Iy,Ia,Iu,rBt,ag+dag,signeo);
    [~,~,~,~,~,~,~,~,Ftm,Ipm,bpm,bpDm,lim,qAm,Sm] = evalFxy(L,Iy,Ia,Iu,rBt,ag-dag,signeo);
    Sg(:,ii) = (Sp - Sm)/(2*dd);
  else
    [~,~,~,~,~,~,~,~,Ftp,Ipp,bpp,bpDp,lip,qAp] = evalFxy(L,Iy,Ia,Iu,rBt,ag+dag);
    [~,~,~,~,~,~,~,~,Ftm,Ipm,bpm,bpDm,lim,qAm] = evalFxy(L,Iy,Ia,Iu,rBt,ag-dag);
  end
  Ftg(:,ii) = (Ftp - Ftm)/(2*dd);
  Ipg(:,ii) = (Ipp - Ipm)/(2*dd);
  bpg(:,ii) = (bpp - bpm)/(2*dd);
  bpDg(:,ii) = (bpDp - bpDm)/(2*dd);
  lig(:,ii) = (lip - lim)/(2*dd);
  qAg(:,ii) = (qAp - qAm)/(2*dd);
end

    
end

function [Fx,FA,FB,rA,zA,Fn,Brn,Bzn,Ft,Ip,bp,bpD,li,qA,S] = evalFxy(L,Iy,Ia,Iu,rBt,ag,signeo)
%% Solve GS to get new flux map
sig = nargout > 14;
Fx = meqFx(L,Iy,[Ia;Iu]);
sIp = sign(sum(Iy(:)));
[rA,zA,FA,dr2FA,dz2FA,drzFA,~,~,~,~,~,~,~,~,...
  ~,~,FB,lB,~,Opy] = meqpdom(Fx,sIp,L.P.isaddl,L);
assert(lB,'could not find boundary');

if ~isempty(L.P.infct)
  [Fn,Brn,Bzn,~,~,~,~] = L.P.infct(L.rx,L.zx,Fx,L.P.rn,L.P.zn,L.inM);
else
  Fn = [];  Brn = []; Bzn = [];
end
  
nB = numel(FB);
sIp = sign(sum(Iy(:)));
[F0,F1] = meqpdomlim(FA,FB,L.nD); % flux limits of various domains

% Total integrals
[~,TpDg,ITpDg] = L.bfct(1,L.bfp,Fx,F0,F1,Opy,L.ry,L.iry);
[Ip,~,~,~,~,~,~,Ft,~,bp,~,~,li,~] = ...
  meqint(L.fPg,L.fTg,TpDg,ITpDg,ag,Fx,Opy,L,rBt);
for iD=1:L.nD
  [~,~,~,~,~,~,~,...
    ~,~,bpD(iD),~,~,~,~] = ...
    meqint(L.fPg,L.fTg,TpDg(iD,:),ITpDg(iD,:),ag,Fx,Opy==iD,L,rBt);
end

% q axis using small diamagnetic approximation (LIUQE 2015 paper eq. 88)
nA = int32(numel(FA));
qA = zeros(nA,1);
[gAg,IgAg] = L.bfct(5,L.bfp,[],F0,F1); IgAg = 2e-7*L.idsx*IgAg;
for iA=1:nA
  TyAg = gAg(iA,:).*(rA(iA)*L.fPg+L.fTg/rA(iA)).'; % TyAg: Axis current (jphi.dR.dZ) contribution from each basis function
  cqA = 4e-7*pi*rA(iA)^2*rBt*sqrt(dr2FA(iA)*dz2FA(iA)-drzFA(iA)^2)/(abs(dr2FA(iA)+dz2FA(iA))*L.dsx);
  qA(iA) = (rBt.^2+(IgAg(iA,:).*L.fTg.')*ag)/(cqA*TyAg*ag);
end


if sig
  [~,~,~,~,FtQ,Ft0Q] = ...
    meqintQ(L,F0,F1,rBt,ag,Fx,Opy,sIp);
  FtPQ  = FtQ+Ft0Q; % total Plasma toroidal flux on pQ grid
  
  % CDE conductivity integral
  S = zeros(1,nB);
  for iD=1:nB
    [~,~,~,TQ] = ...
      meqprof(L.fPg,L.fTg,ag,L.pQ(:).^2,F0,F1,rBt,L.bfct,L.bfp,L.idsx);
    iTsQ = signeo(:,iD)./TQ(:,iD).^2; % sigma/F^2 profile appearing in various places
    S(:,iD) =  rBt/2/pi * trapz(FtPQ(:,iD), iTsQ);
  end
end

end




