function [L,LX,LY] = fge(tok,shot,t,varargin)
%FGE - FORWARD EVOLUTIVE GRAD-SHAFRANOV SOLVER
% Solves the evolutive plasma free boundary equilibrium problem
% Evolves coil and vessel current coupled to plasma GS equation
% [L,LX,LY] = FGE(TOK,SHOT,T,'PAR',VALUE,...)
%
% L = FGE('TOK',SHOT,...) returns only the pre-computed quantities as provided by 
% FGEP<TOK>, FGEG<TOK> and FGEC for later use.
% [L,LX] = FGE('TOK',SHOT,...) additionally returns the data in LX
% generated from FGEX<TOK>
%
% For general help,see README.MD and MEQHELP
% For help on LY,  see MEQT and FBTT
% For help on LX,  see FGEX and FGEX<TOK>
% For help on L.P, see MEQP, MEQP<TOK>, FGEP and FGEP<TOK> 
% For help on L.G, see MEQG, MEQG<TOK>, FGEG and FGEG<TOK> 
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

assert(ischar(tok),'tok must be a string');
assert(isnumeric(shot)||ischar(shot)||iscell(shot),'shot must be numeric (for database) or string (for file)');
if nargin>=3, assert(~mod(numel(varargin),2),'varargin must come in ''parameter'',value pairs'); end

% Tokamak specific functions
tok = lower(tok); % tok: 'tcv' 'ids'
fgeptok = str2func(['fgep' tok]);
fgegtok = str2func(['fgeg' tok]);

% Algorithm configuration parameters, geometry and anciliary data
P = fgeptok(shot(1),varargin{:}); % parameters from sources
G = fgegtok(shot(1),P);  % geometry from sources
P = fgep(P);   % Assemble all parameters
G = fgeg(G,P); % Assemble geometry
L = fgec(P,G); % Consolidate parameters

if nargout >= 2
  assert(isnumeric(t) || ~isempty(t),'time must be a numeric and nonempty');
  % Get inputs data for simulation
  LX = fgex(tok,t,L);
  % compute linearization and initial condition
  if isequal(L.P.algoNL,'all-nl')
    L = fgel(L,meqxk(LX,1));
  end
end

if nargout >= 3
  % Time stepper
  LY = fget(L,LX);
end

end
