%LIUP  LIUQE default algorithm configuration parameters
% P = LIUP(P,'PAR',VAL,...) checks and consolidates the structure P with
% default configuration parameters, optionally replacing or adding parameters
% with specified values.
%
%
% .argout        Optional list of fields for LY (output) to be specified as a cell array of strings with list of VAR=EXPR or VAR, e.g. P.argout={'PpQg','zIp=zA*Ip'}
% .Bmerr     [T] Error in B field measurements
% .chimax        Maximum cost function for accepting equilibrium
% .elomin        If kappa>elomin, regularization weight decreases with elongation, else it is equal to wreg.
% .fdz           Scaling factor for dz equation in least-squares problem
% .Fferr    [Wb] Error in flux measurement
% .fetype        Finite Element type for liu initial guess
% .Fterr [T^2/m] Error in (diamagnetic flux) x (rBt) measurement
% .Iaerr     [A] Error in active coil current measurements
% .idml          If idml~0, use diamagnetic flux measurement
% .Iperr     [A] Error in plasma current measurement
% .ipm           If ipm=1, use Interior-point Method for least-squares fit of currents
% .ipmi          [inactive] Maximum number of iterations in inner loop of IPM for Jh fit.
% .Ipmin     [A] Minimum plasma current to try a reconstruction
% .ipmo          [inactive] Maximum number of iterations in outer loop of IPM for Jh fit.
% .itera         Maximum number of iterations
% .iterfrz       Number of Picard iterations before which plasma domain is frozen to ease convergence (inactive if itert>0)
% .iterh         Maximum number of FE solver iterations
% .iters         Number of block-solver iterations. 0 uses non-iterative solver
% .itert         Maximum number of real-time iterations (0=off-line with itera)
% .Iuerr     [A] Error in vessel current measurements
% .kalman        Kalman filter
% .nelem         Number of FE for initial guess of the current distribution
% .psichco       Convergence criterion on flux
% .rip       [m] Inner boundary for FE
% .rop       [m] Outer boundary for FE
% .robust        If true, fit measurements using robust fitting (requires Statistics and Machine Learning toolbox)
% .slx           If true, liu() calls Simulink version liutsim.m
% .stab          Proportional gain to stabilize vertical position
% .stabcint      Proportional gain to stabilize vertical position. This should be used only if zfct=liuz0.
% .stabder       Derivative gain to stabilize vertical position. This should be used only if zfct=liuz0.
% .stabint       Integral gain to stabilize vertical position
% .stabz         Shift equilibrium to stabilize vertical position
% .tolh          Tolerance for FE solver iterations
% .wBm           Relative weight of B-field measurements
% .wFf           Relative weight of flux measurements
% .wIa           Relative weight of active coil current measurements
% .wIp           Relative weight of plasma current measurements
% .wIu           Relative weight of vessel current measurements
% .wreg          Weight(s) for current regularization terms (scalar, or vector for each entry)
% .wregadapt     Decay rate for regularisation weight with elongation (0 keeps regularisation constant). 
% .zfct          Vertical position stabilisation function handle.
% .z0            Provides user-defined estimate of vertical position. This should be used only if zfct=liuz0.
% .zlp       [m] Lower boundary for FE
% .zup       [m] Upper boundary for FE
%
% .idens         [inactive] if idens=1, LIUQE looks for line integrated density measurements. This should only be used when ifara=1. 
% .ifara         [inactive] if ifara=1, LIUQE looks for Faraday rotation measurements. This should only be used when idens=1. 
% .iinv          [inactive] if iinv>= 2 --> q(0)=1.+float(iinv)/100, if iinv = 1 --> Liuqe reads r,z coordinates of sawtooth inversion points, if iinv = 0 --> q(0) is left floating, If iinv < 0 --> q(0)=1.+float(iinv)/100
% .ipolarimetry  [inactive]
% .ipres         [inactive] if ipres=1, LIUQE looks for pressure measurements. In this case the LIUQE times coincide with the Thomson times. 
% .Ipmeas        [inactive] if Ipmeas=1., LIUQE looks for an independent plasma current measurement. 
%
% Parameters common to MEQ added by meqp.
%
% List of shared parameters with specific default values
%  * gsxe=1 is the default value for LIU
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function P = liup(P,varargin)

 %% Default parameters
 D = meqp();
 % execution parameters
 D.bfct         = @bfabmex;
 D.bfp          = [1 2];
 D.zfct         = @(P,rst,z0,zA,dz) dz;
 D.fdz          = 1e-5;
 D.itert        = 0;
 D.iters        = 48;
 D.iterh        = 0;
 D.tolh         = 0;
 D.dasm         = 2;
 D.xdoma        = 0;
 % configuration
 D.Bmerr        = Inf;
 D.Iaerr        = Inf;
 D.Iuerr        = Inf;
 D.chimax       = 2;
 D.elomin       = 1.2;
 D.Fferr        = Inf;
 D.fetype       = 'bil'; % bilinear finite elements for initial guess
 D.idml         = 1;
 D.ipm          = 0;
 D.ipmi         = 0;
 D.ipmo         = 0;
 D.itera        = 100;
 D.ivesm        = 0; 
 D.itera        = 100;
 D.iterfrz      = 50;
 D.kalman       = false;
 D.nelem        = 6;
 D.Ipmeas       = 1;
 D.Iperr        = Inf;
 D.Ipmin        = 0;
 D.psichco      = 1e-4;
 D.rip          = NaN;
 D.rop          = NaN;
 D.robust       = false;
 D.algoNL       = 'picard';
 D.slx          = 0;
 D.stab         = 1;
 D.stabint      = 1;
 D.stabz        = 1;
 D.stabcint     = 1;
 D.stabder      = 1;
 D.Fterr        = Inf;
 D.wIp          = 1;
 D.wIa          = [];
 D.wIu          = [];
 D.wBm          = [];
 D.wreg         = 0;
 D.wregadapt    = 0;
 D.wFf          = [];
 %D.z0 - Do not specify to keep default behavior
 D.zlp          = NaN;
 D.zup          = NaN;
 
 D.agerr        = 1;
 D.dzerr        = 1;
 D.wag          = 0;
 D.wdz          = 0;
 
 % MEQP parameters with different default value
 D.gsxe         = 1; % Include external current inside the grid in Iy
 
 %% Overwrites extra parameters
 for k = 1:2:length(varargin), P.(varargin{k}) = varargin{k+1}; end
 %% Set defaults for .slx 
 if isfield(P, 'slx') && P.slx, D.bfct= @bf3pmex; D.bfp=false; end
 if isfield(P, 'slx') && P.slx, D.itert= 1; end
 %% Check robust fitting availability
 if isfield(P, 'robust') && P.robust && isempty(which('robustfit'))
  error('MEQP:noRobustFit','P.robust is true but robustfit is unavailable.\nThis option requires the Statistics and Machine Learning toolbox')
 end
 %% Set defaults to call liut as iterator for Newton method
 if isfield(P, 'algoNL') && strcmp(P.algoNL, 'jfnk')
  D.itert = 3;
  D.iters = 0;
 end

 %% Set defaults
 for k = fieldnames(D)'
  if ~isfield(P,k{1}), P.(k{1}) = D.(k{1}); end
 end
 
 %% Check configuration parameters
 msgid = 'LIUQE:ParameterNotSupported';
 msge = 'LIUQE Parameter error: %s %g';
 assert( ~mod(log2(P.nz),1)        ,msgid,msge,'NZ not a power of 2')
 assert( P.Bmerr  >  0             ,msgid,msge,'BMERR=0'            )
 assert( P.Iaerr  >  0             ,msgid,msge,'IAERR=0'            )
 assert( P.Iuerr  >  0 || ~P.ivesm ,msgid,msge,'IUERR=0'            )
 assert( P.Fferr  >  0             ,msgid,msge,'FFERR=0'            )
 assert( P.Iperr  >  0             ,msgid,msge,'IPERR=0'            )
 assert( P.Fterr  >  0             ,msgid,msge,'FTERR=0'            )
 assert( P.agerr  >  0             ,msgid,msge,'AGERR=0'            )
 assert( P.dzerr  >  0             ,msgid,msge,'DZERR=0'            )
 
 assert( all(P.iqR  >   0)         ,msgid,msge,'IQR<=0'             )

 assert(~P.slx || P.itert>0,msgid,msge,...
  'ITERT must be >0 for Simulink case');
 assert(~P.slx || P.iters>0,msgid,msge,...
  'ITERS must be >0 for Simulink case');
 assert(~P.slx || P.ilim==1,msgid,msge,...
  'ILIM must be 1 for Simulink case');
 assert(~P.slx || contains(func2str(P.bfct),{'bf3pmex','bf3imex','bfefmex'}),msgid,msge,...
   'only bf3pmex/bf3imex/bfefmex allowed for bfct if running Simulink case')
 assert(~P.slx || P.gsxe==1,msgid,msge,...
  'GSXE must be 1 for Simulink case');
 k = {'tokamak';'shot'};
 P = orderfields(P,[k;setdiff(fieldnames(P),k,'sorted')]);

end
