function [dLY_dIa,dLY_dIu,dLY_dagcon] = fgslin_full(L,LY)
% [dLY_dIa,dLY_dIa,dLY_dagcon] = fgslin_full(L,LY);
% Linearized static equilibrium.
% Computes dLY/dIa (derivative of LY w.r.t. a coil current)
% and dLY_dagcon (derivative of LY w.r.t. ag constraint quantities)
% Based on finite differences. Counterpart of fgslin for checking etc.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%% linearize around coil currents Ia
LX0 = LY;
dIa = 1e-4*abs(LY.Ip); % Ia perturbations to consider [A]
linfields = fieldnames(LY)'; % fieldnames around which to compute linerization

dLY_dIa = struct();
dLY_dIu = struct();
dLY_dagcon = struct();

L.P.debug=0; % silence these repetitive calls

disp('perturbing equilibria');
% coils
for ia = 1:L.G.na
  fprintf('\nCoil:ia=%d/%d\n',ia,L.G.na)
  LX = LX0;
  
  LX.Ia(ia) = LX0.Ia(ia) + dIa;
  [LYp] = fgst(L,LX);
  LX.Ia(ia) = LX0.Ia(ia) - dIa;
  [LYm] = fgst(L,LX);
  
  % save fields
  dLY_dIa = storederivfield(dLY_dIa,LYp,LYm,linfields,dIa,ia);
end

% vessel
dIu = 1e-4*abs(LY.Ip); % Iu perturbations to consider [A]
for iu = 1:L.G.nu
  fprintf('\nVessel:iu=%d/%d\n',iu,L.G.nu)
  LX = LX0;
  
  LX.Iu(iu) = LX0.Iu(iu) + dIu;
  [LYp] = fgst(L,LX);
  LX.Iu(iu) = LX0.Iu(iu) - dIu;
  [LYm] = fgst(L,LX);
  
  % save fields
  dLY_dIu = storederivfield(dLY_dIu,LYp,LYm,linfields,dIu,iu);
end

% Derivatives w.r.t. quantities in LX that define ag constraints (via meqagcon)
for iagcon = 1:numel(L.P.agcon)
  myagcon = L.P.agcon{iagcon}; % ag constraining quantity to perturb
  fprintf('\nagcon:%s\n',myagcon)
  LX = LX0;
  dagcon = LX0.(myagcon)*1e-3;
  LX.(myagcon) = LX0.(myagcon) + dagcon;
  [LYp] = fgst(L,LX);
  
  LX.(myagcon) = LX0.(myagcon) - dagcon;
  [LYm] = fgst(L,LX);
  
  % save fields
  dLY_dagcon = storederivfield(dLY_dagcon,LYp,LYm,linfields,dagcon,iagcon);
  
end
end


function dLY_dX = storederivfield(dLY_dX,LYp,LYm,fields,h,ia)
for field = fields
  myfield=field{:};
  dfield_dIa = (LYp.(myfield)-LYm.(myfield))/(2*h);
  if isscalar(LYp.(myfield)) || iscolumn(LYp.(myfield))
    % scalar or vector
    dLY_dX.(myfield)(:,ia) = dfield_dIa;
  elseif ismatrix(LYp.(myfield))
    dLY_dX.(myfield)(:,:,ia) = dfield_dIa;
  elseif ndims(LYp.(myfield)) == 3
    dLY_dX.(myfield)(:,:,:,ia) = dfield_dIa;
  else
    error('invalid field dimensionality for field %s',myfield)
  end
end
end
