%FBTPTCV  TCV FBT configuration parameters
% P = FBTPTCV(SHOT,'PAR',VAL,...) returns a structure P with configuration
% parameters from the PCS tree for SHOT, optionally replacing or adding
% parameters with specified values. See also FBTP.
% Most parameters are inherited from MGAMS. For some documentaion see:
%  https://spcwiki.epfl.ch/wiki/Mgams_help
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function P = fbtptcv(shot,varargin)
 P = meqptcv(shot);
 
 assert(~isempty(which('mgp')),'mgp.m not found - can''t load shot MGAMS programmed data for FBT-TCV');
 % Load from PCS tree
 mdsconnect(P.mdsserver); % make sure to connect to tcvdata because mgp does not necessarily do it
 Pmgp = mgp(shot);

 % Shots prior to 56469 did not have gp* variables and these will not have the correct number of columns
 if (shot > 0 && shot < 56469)
   Pmgp.fix;
 end

 for ifield=(Pmgp.NAME).'
   if startsWith('debug_p',ifield{1}), continue; end
   P.(ifield{1}) = Pmgp.(ifield{1});
 end
 
 % Process additional arguments from debug_p1
 if ~isprop(Pmgp,'debug_p1')
  warning('FBTE:fbtptcv:debug_p1_missing','debug_p1 is missing from mgp object, check the version of shotdesign used');
 else
  P1 = Pmgp.debug_struct(1);
  for ifield = fieldnames(P1).'
   P.(ifield{1}) = P1.(ifield{1});
  end
 end
 
 % kA to A, Wb/rad to Wb and some defaults
 P.shot    = mdsdata('$SHOT');
 P.placu1  = 1000*P.placu1;
 P.dissi   = 1e-6*P.dissi ;
 P.strki   = 1e-6*P.strki ;
 P.t       = 0.01+P.timefac*(P.timeeq-P.timeeq(1));
 P.bfct    = @bffbt;    % FBT basis function set
 P.sela    = {'E','F'}; % No OH or G coils in FBTE
 P.selu    = 'n';       % no vessel in FBTE
 P.selx    = 'XF';      % Historical FBTE grid
 P.zl      = -P.zu;     % New parameter for asymmetric grids
 
 P.agfitfct = @fbtfit;  % TCV legacy fbt fit function

 % adds/overwrites parameters
 for k = 1:2:nargin-1
  P.(varargin{k}) = varargin{k+1};
 end
 
 % Cannot use a different static tree when preparing the MODEL (-1) shot
 assert(shot ~= -1 || P.static == -1,'Cannot specify a custom static tree when preparing the MODEL (-1) shot');
 
 if ~all(isfield(P,{'gpvrr','gpvrz','gpvzz','gpve','gpvd'}))
  if any(isfield(P,{'gpvrr','gpvrz','gpvzz','gpve','gpvd'}))
   msg = {'Some of the gpv* variables are missing, replacing all with NaNs.',...
          'Please provide ''gpvrr'', ''gpvrz'', ''gpvzz'', ''gpve'' and ''gpvd''.'};
   warning('fbtptcv:gpv_incomplete',strjoin(msg,'\n'));
  end
  P.gpvrr = NaN(size(P.gpr));
  P.gpvrz = NaN(size(P.gpr));
  P.gpvzz = NaN(size(P.gpr));
  P.gpve  = NaN(size(P.gpr));
  P.gpvd  = NaN(size(P.gpbd));
 end
 
 % Treat Ip and B0 signs
 % NOTE: We assume here that in MGAMS every GP variable is set for Ip,B0>0
 assert(all(P.placu1>0) && all(P.bzero>0),'placu1,bzero should be >0 from MGAMS')
 P.placu1  = P.iohfb*P.placu1;
 P.gpia    = P.iohfb*P.gpia;
 P.gpfa    = P.iohfb*P.gpfa;
 P.gpbr    = P.iohfb*P.gpbr;
 P.gpbz    = P.iohfb*P.gpbz;
 P.gpba    = (1+sign(P.iohfb))*pi+P.gpba;
 P.gpcr    = P.iohfb*P.gpcr;
 P.gpcz    = P.iohfb*P.gpcz;
 P.gpvrr   = P.iohfb*P.gpvrr;
 P.gpvrz   = P.iohfb*P.gpvrz;
 P.gpvzz   = P.iohfb*P.gpvzz;
 P.bzero   = P.if36fb*P.bzero;
 P.qzero   = P.iohfb*P.if36fb*P.qzero;
 
 % Use unified parameters
 P = fbtptcv_update(P);
 
 % Check consistency of sizes between pre-existing GPs
 P = gpcheck(P);
 
 % Convert old shape parameters
 n = max(P.ilia(logical(P.iansha)));
 P.rlia1(end+1:n,:) = NaN;
 P.zlia1(end+1:n,:) = NaN;
 for k = find(P.iansha)
  n = P.ilia(k);
  x = 2*pi/n*(0:n-1)';
  P.rlia1(:,k) = NaN;
  P.zlia1(:,k) = NaN;
  P.rlia1(1:n,k) = (P.rmajo1(k)+P.rmino1(k)*cos(x+P.delta1(k)*sin(x)-P.hlamd1(k)*sin(2*x)));
  P.zlia1(1:n,k) = (P.zmajo1(k)+P.rmino1(k)*P.cappa1(k)*sin(x));
 end
 %   fbtgp(P,r      ,z      ,b,fa,fb,fe,br,bz,ba,be,cr,cz,ca,ce,vrr,vrz,vzz,ve)
 P = fbtgp(P,P.rlim1,P.zlim1,1, 0, 1, 0,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
 P = fbtgp(P,P.rlia1,P.zlia1,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
 P = fbtgp(P,P.rbro ,P.zbro ,0,[],[],[], 0,[],[], 0,[],[],[],[],[] ,[] ,[] ,[]);
 P = fbtgp(P,P.rbzo ,P.zbzo ,0,[],[],[],[], 0,[], 0,[],[],[],[],[] ,[] ,[] ,[]);
 % Adjust selection for plasma initial guess
 P.gpb = gpb(P);
 
 % Coil names
 dima = [cellstr(num2str((1:8).','E_%03d'));cellstr(num2str((1:8).','F_%03d'))];
 
 % Active PFC currents constraints
 mask = isfinite(P.gpie.*P.gpid);
 
 % Avoid conflicts between icoilon and gpia
 conflicts = mask & ~P.icoilon;
 if any(conflicts(:))
   % Prepare error message with lists of times (runs) and coils with conflicts
   runs_c = find(any(conflicts,1));
   msg = cell(1,numel(runs_c)+1);
   msg{1} = 'Conflicts between icoilon and gpia settings.';
   for irun = 1:numel(runs_c)
     msg{1+irun} = sprintf('  %d: %s',irun,strjoin(dima(conflicts(:,runs_c(irun)))));
   end
   error('fbtptcv:IaConflict',strjoin(msg,'\n'));
 end

 % translate isaddl to specific doublet parameters
 if ~isfield(P,'idoublet')
   P.idoublet = deal(~logical(P.isaddl));
 end
 
 % Combine with old dissipation of PFC currents
 P.gpia(~mask) = 0;                                           % if no active constraint, set value to 0
 P.gpie(~mask) = 1;                                           % Default weight
 P.gpie(1:8,:) = P.gpie(1:8,:)*sqrt(2)*18/17;                 % Increase default weight for E coils
 P.gpie(~mask & ~P.icoilon) = P.gpie(~mask & ~P.icoilon)/1e3; % Coils "deselected" have large weights for their current.
 
 % Convert old weights for PFC currents dipole
 % Wij = Wpsi*strki/dij^2*(w(i)*Ia(i) - w(j)*Ia(j))^2
 % Variables: dij -> dc, w(i) -> gpdw
 dima = [cellstr(num2str((1:8).','E_%03d'));cellstr(num2str((1:8).','F_%03d'))];
 meqmdsopen(P.static,'STATIC',[],P.mdsserver);
 rc = mdsdata('STATIC("R_C")[$1]',dima);
 zc = mdsdata('STATIC("Z_C")[$1]',dima);
 wc = mdsdata('STATIC("W_C")[$1]',dima);
 hc = mdsdata('STATIC("H_C")[$1]',dima);
 % E coils have 4 vert. stacked filaments
 % F coils have 2 rad.  stacked filaments
 Ec=1:8;Fc=9:16;
 rff = [rc(Ec);rc(Fc)-wc(Fc)/4];   % R of first coil filament
 rfl = [rc(Ec);rc(Fc)+wc(Fc)/4];   % R of last  coil filament
 zff = [zc(Ec)-hc(Ec)*3/8;zc(Fc)]; % Z of first coil filament
 zfl = [zc(Ec)+hc(Ec)*3/8;zc(Fc)]; % Z of last  coil filament
 dc = sqrt((rff(2:end)-rfl(1:end-1)).^2 + (zff(2:end)-zfl(1:end-1)).^2); % Distance from the last filament to the first of the next coil
 nt = numel(P.t);
 dw = [17/2*ones(8,1);18*ones(8,1)]; % Weights for each coil current (from current in each turn to current in each filament)
 P.gpdw = repmat(dw,1,nt);
 P.gpde = repmat(dc,1,nt); % Weight is inversely proportional to distance.
 
 % Cost function weight
 [P.gpfd,P.gpbd,P.gpcd] = deal(repmat(5e-3,1,P.nruns));
 P.gpid = P.gpfd./(36*pi*sqrt(2*P.dissi));
 P.gpdd = P.gpfd./( 2*pi*sqrt(  P.dipol));
 
 % Coil limits
 if shot>0 && shot<1e5, base_shot = shot;
 elseif shot>1e6,       base_shot = floor(shot/1e3);
 else,                  base_shot = -1;
 end
 meqmdsopen(base_shot,'BASE',[],P.mdsserver);
 P.limc = [eye(16);mdsdata('\COIL_PROT:CEQ[1..8,[$1]]',dima)];
 P.liml = NaN(size(P.limc,1),1); P.limu = P.liml;
 P.liml(17:end) = mdsdata('\COIL_PROT:ILOW[1..8]');
 P.limu(17:end) = mdsdata('\COIL_PROT:IHIG[1..8]');
 meqmdsopen(P.static,'STATIC',[],P.mdsserver);
 P.limu(1:16) = mdsdata('STATIC("INOM_A")[$1]',dima);
 P.liml(1:16) = -P.limu(1:16);
 meqmdsopen(shot,'PCS',[],P.mdsserver);
 
 % Grid choice
 validatestring(P.selx,{'','X','XF'},mfilename,'P.selx');
end

function P = fbtptcv_update(P)
 % update TCV parameters into FBT/MEQ ones
 
 % quantities which come from mds tree for legacy reasons but
 % should be loaded in fbtxtcv
 P.Ip  = P.placu1;    
 P.rBt = P.bzero*0.88;
 if ~isfield(P,'qA')
   P.qA  = P.qzero;
 end
 
 % quantities which may be provided via fbt(...) calls
 % For these, if no user-provided values are given, write TCV tree values
 if ~isfield(P,'dipol'), P.dipol = P.strki;  end
 if ~isfield(P,'niter'), P.niter = P.itamax; end
 if ~isfield(P,'tol'),   P.tol   = P.testa;  end

 % Remove legacy field names for clarity
 P = rmfield(P,{'qzero','bzero','placu1','strki','itamax','testa'});
end

function P = gpcheck(P)

 nruns = P.nruns;

 % Since the mgp class removes final rows with only NaNs in time-dependent
 % arrays, we need here to restore these so that each GP parameter has the
 % same number of rows.

 list = strcat('gp',{'r','z','b','fa','fb','fe','br','bz','ba','be','cr','cz','ca','ce'});
 nvar = numel(list);
 
 % Get number of rows for each variable and check number of columns matches nruns
 rows = NaN(nvar,1);
 for ii = 1:nvar
  sz = size(P.(list{ii}));
  assert(sz(2) == nruns,'FBTE:Inconsistent number of columns for %s. Expected %d found %d',list{ii},nruns,sz(2));
  rows(ii) = sz(1);
 end
 
 % Now make sure all of them have the same number of rows
 nrows = max(rows);
 for ii = 1:nvar
  P.(list{ii})(rows(ii)+1:nrows,:) = NaN(nrows - rows(ii),nruns);
 end
 
 % Check that we've made things right
 rows0 = size(P.(list{1}),1);
 for ii = 2:nvar
  assert(size(P.(list{ii}),1) == rows0, 'FBTE:Inconsistent number of rows between gp%s and gp%s',list{1},list{ii});
 end
 
 % Check dimensions
 assert(size(P.gpia,1)==16,'FBTE:gpia must have 16 rows')
end

function gpb = gpb(P)
 % Attempt to eliminate points defining strike points characteristics
 % from initial plasma guess (typically for old TCV shots). 
 gpb = P.gpb;
 for k = 1:P.nruns
  % X-points
  ibro = P.gpbr(:,k)==0;
  ibzo = P.gpbz(:,k)==0;
  rzbro = [P.gpr(ibro,k),P.gpz(ibro,k)];
  rzbzo = [P.gpr(ibzo,k),P.gpz(ibzo,k)];
  iX = ismember(rzbro,rzbzo,'rows');
  rX = rzbro(iX,1);
  zX = rzbro(iX,2);
  % Initiate loop
  bo = gpb(:,k) == 1 & ~isnan(P.gpr(:,k)) & ~isnan(P.gpz(:,k));
  kk = 0;
  while kk < 10
   % Axis as center of bounding box of points with gpb=1
   bb = [min(P.gpr(bo,k)) max(P.gpr(bo,k)) min(P.gpz(bo,k)) max(P.gpr(bo,k))]*0.5;
   rA1 = bb(1) + bb(2); zA1 = bb(3) + bb(4);
   % Axis as barycenter of active control points
   no = sum(bo);
   rA2 = sum(P.gpr(bo,k))/no; zA2=sum(P.gpz(bo,k))/no;
   % Weight both answers (1 for bounding box, 3 for barycenter)
   rA = (rA1 + 3*rA2)/4; zA = (zA1 + 3*zA2)/4;
   % Set gpb=0 for points outside of so-defined X-point polygon
   b = bo & bavxmex(rX,zX,rA-rX,zA-zX,0,P.gpr(:,k),P.gpz(:,k));
   kk = kk+1;
   % If gpb did not change, stop
   if isequal(b,bo), break;end
   bo = b;
  end
  gpb(:,k) = b;
 end
end
