function [Ei,Vi,LYi] = fgeeig(L,varargin)
% FGEEIG Calculate unstable eigenvalues of system using fgess
% 
% [Ei,Vi,LYi] = fgeeig(L,varargin)
% Returns unstable eigenvalues (Ei) of system returned by fgess(L)
% eigenmodes (Vi) and array of perturbed equilibrium information (LYi)
% 
% optional parameters-value pairs: 
%  nmodes: number of modes to compute, default: most unstable
%  imodes: specific mode indices to plot
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

assert(isfield(L,'Mee'),'missing Mee field in L, is this an fge structure?')
assert(isfield(L,'lin'),'missing lin field in L, did you run fgel()?')

p=inputParser;
p.addParameter('nmodes',1); % number of modes to plot
p.addParameter('imodes',[]); % number of modes to plot
parse(p,varargin{:}); P = p.Results;

if isempty(P.imodes)
  nmodes = P.nmodes; imodes = 1:nmodes; % plot some modes
else
  imodes = P.imodes; nmodes = numel(imodes); % plot specific modes
end

% Compute state-space system corresponding to linearized equilibirum stored
% in L
sys = fgess(L);

[V,D] = eig(sys.A);
[E,isort] = esort(diag(D));

%%
% init
Ei = zeros(1,nmodes);
Vi = zeros(size(V,1),nmodes);

for jj=1:nmodes
  ieig = imodes(jj);
  Ei(jj) = E(ieig);
  if L.P.debug>0
    fprintf('growth rate %2.2f\n',E(ieig));
  end
  % Find dIy response to eigenvector direction  
  if L.np==0
    dIydV = L.lin.dIydIe; % state is just coils
  else
    % Se are solving a CDE as well, so state also contains Ip
    iagconIp = contains(L.agconc(:,3),'Ip'); % find index of constraint with Ip
    dIydV = [L.lin.dIydIe,L.lin.dIydCo(:,iagconIp)]; % add Ip contribution
  end
  Veig = V(:,isort(ieig)); % eigenvector
  % Adjust eigenvector sign to make it unambiguous
  Veig = Veig*sign(sum(L.zzy(:).*(dIydV*Veig)));
  Vi(:,jj) = Veig; % store

  dIy = reshape(dIydV*Veig,L.nzy,L.nry);
  
  if nmodes>1
    subplot(1,nmodes,jj); cla;
  end
  
  dIe = Veig(1:L.ne);
  if isfield(L.G,'Mza')
    LY.Fz = reshape([L.G.Mza,L.G.Mzu]*dIe + L.G.Mzy*dIy(:),L.nzz,L.nrz);
  end
  LY.Fx = meqFx(L,dIy,dIe);
  LY.Ia = dIe(1:L.G.na); 
  LY.lB = false; [LY.rX,LY.zX] = deal([]);
  LY.Iy = reshape(dIy,L.nzy,L.nry);
  LY.Iu = dIe(L.G.na+(1:L.G.nu)); LY.Iv = L.G.Tvu*dIe(L.G.na+(1:L.G.nu));
  LY.Iy = dIy;
  LY.Ip = 0; % stresses that this is not a real plasma
  LY.t  = NaN;
  LYi(jj) = LY; %#ok<AGROW>
end
