%LIUPCREATE  CREATE(ITER) LIUQE algorithm configuration parameters
% P = LIUPCREATE(EQUIL,'PAR',VAL,...) returns a structure P with algorithm
% configuration parameters for EQUIL, optionally replacing or adding parameters
% with specified values. See also LIUP.
% .selfm        Regular expression for sensor selection
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function P = liupcreate(equil,varargin)

%% Load equilibrium file
file = createf(equil);
if isempty(file)
  error('Cannot find file for equilibrium #%05d',equil)
else
  GeomFile = file{1};
end

if numel(varargin)==0, P = struct; end % init

% Default parameters, not defined by LIUQE xml file
P.tokamak   = 'LIUQE-ITER(C)';
P.tok       = 'create';
P.shot      = equil;
P.fetype    = 'tri'; % triangular basis function elements
P.selu      = 'n'; % default vessel representation: none
P.vesselmode = 0;  % default vessel treatment mode

% Get external parameters if any
for k = 1:2:numel(varargin)
  P.(varargin{k}) = varargin{k+1};
end

if ~isfield(P,'liupxml') || isempty(P.liupxml)
  fprintf('no xml source selected, using default liupcreate parameters\n')
  P = liupcreate_m(P,GeomFile);
else
  P = liupcreate_XML(P,GeomFile);
end

% Override or add additional parameters (if needed)
for k = 1:2:numel(varargin)
  P.(varargin{k}) = varargin{k+1};
end

end

function P = liupcreate_XML(P,GeomFile)
% Load liuqe parameters for CREATE from XML files
xmlLIUQEfile  = P.liupxml;   % LIUQE-specific configuration file
xmlConfigFile = P.confxml;  % General equilibrium reconstruction configuration

assert(~isempty(which(xmlLIUQEfile)), 'could not find file %s',xmlLIUQEfile );
assert(~isempty(which(xmlConfigFile)),'could not find file %s',xmlConfigFile);
assert(~isempty('XML2meqp'),'could not find file XML2meqp(). Please add CREATE OPE883 interface scripts to path');
PXML = XML2meqp(xmlLIUQEfile); % get LIUQE code parameters
% append to previous struct
for field=fieldnames(PXML)'
  P.(field{:}) = PXML.(field{:});
end
D = load(GeomFile); % Get default parameters from Geometry file (all sensors)

% Update sensor weight parameters based on those selected in xml configuration file
XMLConfig = xml_read(xmlConfigFile);

%%%% Fix DML flux for retro-compatibility
if not(isfield(D,'diamagneticflux'))
  D.diamagneticflux.names = [];
elseif not(isfield(D.diamagneticflux,'names'))
  dval = D.diamagneticflux;
  D = rmfield(D,'diamagneticflux');
  [~,~,D.diamagneticflux] = XML2sensNames('ITERmagneticDiagnostics.xml',XMLConfig); % This should be changed so as to use a general database and not an hard-coded one
  D.diamagneticflux.names = D.diamagneticflux.names(1);
  D.diamagneticflux.values = dval;
end
%%%%

[fs,sl,dmloop] = XML2sensWeights([],XMLConfig,D.field_sensors,D.saddle,D.diamagneticflux); %#ok<ASGLU>
P.wBm(:) = fs.weights';
P.wFf(:) = [sl.weight]';
P.wIp(:) = XMLConfig.InputPlasmaCurrentWeight;

% At this point, G.dima is not known yet - therefore we hard-code the final
% coil ordering here to extract the coil current weights in the correct order.
cnames = {'CS3U','CS2U','CS1','CS2L','CS3L','PF1','PF2','PF3','PF4','PF5','PF6','VS'};
P.wIa = zeros(numel(cnames),1);
for ifield = 1:numel(cnames)
  P.wIa(ifield) = XMLConfig.InputCircuitWeights.(cnames{ifield});
end

%% Geometry
if not(isfield(D,'R0'))
  P.r0    = 6.384; % default for ITER
else
  P.r0    = D.R0;
end
P.zlp     = 1/10*(max(D.limiter.z)-min(D.limiter.z))+min(D.limiter.z); % Exclude divertor region from first guess


%% Psin values for flux contouring
temp = strfind({XMLConfig.OutputParameters.FluxContours.FluxContour.CONTENT}, 'psin');
psiidx = 1;
while isempty(temp{psiidx}), psiidx=psiidx+1; end
psin = XMLConfig.OutputParameters.FluxContours.FluxContour(psiidx).ATTRIBUTE;
if psin.compute
  P.pq = sqrt(str2num(psin.position)); %#ok<ST2NM>
else
  P.pq = [];
end

%% Rn values for flux contouring
temp = strfind({XMLConfig.OutputParameters.FluxContours.FluxContour.CONTENT}, 'rn');
ridx = 1;
while isempty(temp{ridx}), ridx=ridx+1; end
ras = XMLConfig.OutputParameters.FluxContours.FluxContour(ridx).ATTRIBUTE;
if ras.compute
  P.raS = str2num(ras.position); %#ok<ST2NM>
else
  P.raS = [];
end

%% V loop computation
if XMLConfig.OutputParameters.Vloop.ATTRIBUTE.compute && not(isempty(XMLConfig.OutputParameters.Vloop.ATTRIBUTE.r))
  P.rn = str2num(XMLConfig.OutputParameters.Vloop.ATTRIBUTE.r); %#ok<ST2NM>
  P.zn = str2num(XMLConfig.OutputParameters.Vloop.ATTRIBUTE.z); %#ok<ST2NM>
  assert(numel(P.rn)==numel(P.zn),'r and z coordinates for vloop computation should have the same number of elements');
  P.infct=@qintmex; 
end

%% Flags for storing output
P.bpflag   = XMLConfig.OutputParameters.PoloidalBeta.ATTRIBUTE.compute;
P.liflag   = XMLConfig.OutputParameters.InternalInductance.ATTRIBUTE.compute;
P.qflag    = XMLConfig.OutputParameters.qProfile.ATTRIBUTE.compute;
P.qminflag = XMLConfig.OutputParameters.qmin.ATTRIBUTE.compute;
q15flag    = XMLConfig.OutputParameters.q15.ATTRIBUTE.compute;
q2flag     = XMLConfig.OutputParameters.q2.ATTRIBUTE.compute;

assert(~P.qminflag || P.qflag, 'must select computation of q profile if qmin output is desired');
assert(~q15flag    || P.qflag, 'must select computation of q profile if q2/1 output is desired');
assert(~q2flag     || P.qflag, 'must select computation of q profile if q3/2 output is desired');

% select 1/q value of rational surfaces to track
P.iqR = []; % init empty
if q15flag, P.iqR = [P.iqR,2/3]; end
if q2flag,  P.iqR = [P.iqR,1/2]; end

end
  

function P = liupcreate_m(P,GeomFile)
% legacy parameter loading function without XML - only this one works in
% stand-alone with meq suite
D = load(GeomFile);

%% adds/overwrite parameters
P.Bmerr   =    0.003;
P.Fferr   =    0.1  ;
P.Iaerr   =    0.01 ;
P.Iuerr   =    0.01 ;
P.Fterr   =    0.1  ;
P.Iperr   =    1000 ;
P.wIa     =    ones(1,length(D.coils_data.turns)-2);
P.wIu     =       [];
P.wreg    =        0;
if not(isfield(D,'R0'))
  P.r0    =    6.384; % default for ITER
else
  P.r0    =     D.R0;
end
P.nelem   =       10;
P.iterh   =       20;
P.iters   =        0;     % use full solver (does not work for Simulink case)
P.tolh    =     1e-3;
P.zlp     = 1/10*(max(D.limiter.z)-min(D.limiter.z))+min(D.limiter.z); % Exclude divertor region from first guess
P.psichco =     1e-6;
P.nz      =       64;
P.idml    =        0;
P.selfm   = '55\.A[ABDEL].*';
P.npq     =       20; % require 20 contour calculations by default

end
