function LY = fgetkl(L,LX,zLY)
% Step linear system given in state space form
% Return usual LY structure
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%% Get current state
zxt = [zLY.Ia;zLY.Iu];
if L.icde, zxt = [zxt;zLY.Ip]; end

%% Update time varying inputs
% Get constraints current time step
Co = LX2Co(L,LX);
zCo = LX2Co(L,zLY);

% Extract constraints which do not belong to current diffusion
if any(L.icde)
  icoin = find(~contains(L.agconc(:,3), 'Ip'));
  Co = Co(icoin);
  zCo = zCo(icoin);
end

% Pack time variable inputs u = {Va,Co,dCodt}; See fgess.m for more details.
ut = zeros(size(L.lin.sysB,2),1);
if ~L.icde
  ut(L.lin.iut) = [Co;(Co-zCo)/LX.dt;LX.Va];
else
  ut(L.lin.iut) = [Co;(Co-zCo)/LX.dt;LX.Va;LX.Ini];
end

%% Step system
xt = L.lin.sysA*zxt + L.lin.sysB*ut;

%% Compute post processing
yt = L.lin.sysC*xt + L.lin.sysD*ut;

%% Add bias
yt = yt + L.lin.yo;

%% Pack outputs into standard LY structure
LY = zLY; % Initialize from previous time slice

% Extracting indices using L.lin.sys.OutputGroup directly is slower
% It is significantly faster to copy the structure first
for ii = fieldnames(L.lin.OutputGroupString)'
  LY.(ii{:}) = reshape(yt(L.lin.OutputGroupString.(ii{:})),meqsize(L,ii{:}));
end

% Set to NaN all the fields which are not computed for the linear post processing
for field = setdiff(fieldnames(zLY),[fieldnames(L.lin.OutputGroupString);'isconverged'])'
  if islogical(LY.(field{:}))
    LY.(field{:})(:) = false;
  else
    LY.(field{:})(:) = NaN;
  end
end

%% Add extra inputs not available in linearization outputs
LY.t = LX.t;
LY.Va = LX.Va;
LY.rBt = LX.rBt;

%% Add flux surface post processing
if L.P.iterq
  % Replace D quantities for single domain for meqposq
  if L.nD==1
    for ii=fieldnames(LY)'
      str =  ii{:};
      if contains(str,'D')
        LY.(str) = LY.(str(1:end-1));
      end
    end
  end
  [LY.rA,LY.zA,LY.FA,LY.dr2FA,LY.dz2FA,LY.drzFA,LY.rX,LY.zX,LY.FX,~,~,~,~,~,...
    LY.rB,LY.zB,LY.FB,LY.lB,LY.lX,LY.Opy] = meqpdom(LY.Fx,sign(LY.Ip),L.P.isaddl,L);
  [F0,F1] = meqpdomlim(LY.FA,LY.FB,L.nD); 
  [LY.PpQ,LY.TTpQ,LY.PQ,LY.TQ,LY.iTQ,LY.PpQg,LY.TTpQg] = ...
  meqprof(L.fPg,L.fTg,LY.ag,L.pQ(:).^2,F0,F1,LX.rBt,L.bfct,L.bfp,L.idsx,L.smalldia);

  LY.nA=int32(numel(LY.rA)); LY.nB = int32(numel(LY.FB)); LY.nX = int32(numel(LY.FX)); % number of boundaries = number of active domains

  LY.aq = zLY.aq;
  LY.aW = zLY.aW;
  [LY] = meqpostq(L,LY);
  
end

%% TODO make some debugging plot

%% Add sanity check for the state

end