%MEQGTCV  TCV MEQ geometry and electromagnetic parameters
% MEQGTCV(P,A,S,L,F,M,X,C,U) returns a structure with geometry and
% electromagnetic parameters for TCV SHOT. 
%   P contains parameter (e.g. 'shot')
%   A selects the active coil sets amongst {'E' 'F' 'G' 'OH'}. 
%   S specifies the number of vessel segments, 
%   L specifies the number of limiter points, 
%   F, M the indices of the flux loops and  magnetic probes, F=-1 for 3D flux loops, 
%   X is one of {'' 'X' 'XF'} and selects the appropriate mesh. 
%   C xspecifies whether to return 'coil' geometry information (true/false)
%   U specifies type of vessel description to use ('v','s','e','n')
%   
% See also LIUGTCV,FBTGTCV,MEQG.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function G = meqgtcv(P,sela,sels,sell,self,selm,selx,selu,selc,varargin)
 
 % Open static tree
 meqmdsopen(P.static,'STATIC',[],P.mdsserver);
 
 %% Dimension/system selection
 % Default
 if isempty(sela), sela = {'E' 'F' 'G' 'OH'}; end
 if isempty(sels), sels =  38; end, dims = 1:sels;
 if isempty(sell), sell = 128; end, diml = 1:512/sell:512;
 if isempty(self) , dimf = 1:38; self = 'F' ;
 elseif self == -1, dimf = 1:38; self = 'F1';
 else             , dimf = self; self = 'F' ; 
 end
 if isempty(selm) , selm = 1:38; end, dimm = selm;
 if isempty(selc) , selc = false;end
 if isempty(selu) , selu='s'; end
 
 % Dims and sel
 G.dima = {};
 for k = 1:length(sela)
  switch sela{k}
   case 'E',   l = 1:8;
   case 'F',   l = 1:8;
   case 'G',   l = 1;
   case 'OH',  l = 1:2;
   case 'TOR', l = 1;
   otherwise, error('unknown coil %s',sela{k});
  end
  G.dima = [G.dima; cellstr(num2str(l(:),[sela{k} '_%03d']))];
 end
 G.dims = cellstr(num2str(dims','S_%03d'));
 G.dimf = cellstr(num2str(dimf','%03d'  ));
 G.dimm = cellstr(num2str(dimm','%03d'  ));
 G.diml = cellstr(num2str(diml','%03d'  ));
 G.dimv = cellstr(num2str((1:mdsdata('STATIC("N_V")'))','v%03d'  ));
 
 %% Active coils
 Tca = mdsdata('STATIC("T_C_A")[*,$1]',G.dima);
 G.rw = mdsdata('STATIC("R_W")');
 G.zw = mdsdata('STATIC("Z_W")');
 G.ww = mdsdata('STATIC("A_W")'); % Conductor width
 G.hw = mdsdata('STATIC("B_W")'); % Conductor height
 G.dw = mdsdata('STATIC("D_W")'); % Conductor diameter (for circular windings)

 G.Twa = mdsdata('STATIC("T_W_C")')*Tca;
 k = any(G.Twa,2);
 G.rw  = G.rw (k  );
 G.zw  = G.zw (k  );
 G.ww  = G.ww (k  );
 G.hw  = G.hw (k  );
 G.dw  = G.dw (k  );
 G.Twa = G.Twa(k,:);
 
 %% Vessel
 G.rv  = mdsdata('STATIC("R_V")');
 G.zv  = mdsdata('STATIC("Z_V")');
 G.wv  = mdsdata('STATIC("A_V")'); % Conductor width
 G.hv  = mdsdata('STATIC("B_V")'); % Conductor height
 G.nv  = mdsdata('STATIC("N_V")');
 
 %% Limiter contour
 G.rl = mdsdata('STATIC("R_T")[$1]',G.diml);
 G.zl = mdsdata('STATIC("Z_T")[$1]',G.diml);
 
 %% Magnetics
 G.rm = mdsdata('STATIC("R_M"  )[$1]',G.dimm);
 G.zm = mdsdata('STATIC("Z_M"  )[$1]',G.dimm);
 G.am = mdsdata('STATIC("ANG_M")[$1]',G.dimm);
 G.rf = mdsdata('STATIC("R_F"  )[$1]',G.dimf);
 G.zf = mdsdata('STATIC("Z_F"  )[$1]',G.dimf);
 
 %% Magnetics Green's functions
 G.Mfa   = mdsdata('STATIC("MUT_"//$1//"_A")[$2,$3]',self,G.dimf,G.dima);
 G.Bma   = mdsdata('STATIC("BT_M_A" )       [$1,$2]',     G.dimm,G.dima);
 
 %% Mesh from STATIC tree
 if ~isempty(selx)
  G.rx = mdsdata('STATIC("R_"//$1//":AXIS")',selx);
  G.zx = mdsdata('STATIC("Z_"//$1//":AXIS")',selx);
  % Green's functions
  G.Mxa   = mdsdata('STATIC("MUT_"//$1//"_A"   )[ *,$2]',     selx,G.dima);
 
  G.Mfx   = mdsdata('STATIC("MUT_"//$1//"_"//$2)[$3, *]',self,selx,G.dimf);
  G.Bmx   = mdsdata('STATIC("BT_M_"//$1        )[$2, *]',     selx,G.dimm);
 end
 
 %% Coil locations (for PLOTS)
 if selc
  G.dimc = mdsdata('STATIC("DIM_C")');
  G.Tca  = mdsdata('STATIC("T_C_A")[*,$1]',G.dima);
  G.rc   = mdsdata('STATIC("R_C")');
  G.zc   = mdsdata('STATIC("Z_C")');
  G.wc   = mdsdata('STATIC("W_C")');
  G.hc   = mdsdata('STATIC("H_C")');
 end
 
 %% Vessel description
 if isequal(selu,'n')
   dimu=[];
 else
   
   nu = mdsvalue('STATIC("N_"//$1//"")',selu); % default nu for this vessel description
   if isfield(P,'nu') && ~isempty(P.nu)
     % only 'e' allows smaller custom nu
     assert((contains(selu,'e') && P.nu<=nu) || P.nu==nu,'invalid nu %d for this selu %s',nu,selu);
     nu = P.nu; % use custom nu if ok
   end
   
   % dimu for retrieving vessel elements/modes from static tree
   dimu = mdsdata('STATIC("DIM_"//$1)',selu);
   dimu = dimu(1:nu);
   
   if (selu(1) == 'v')
    G.Tvu = eye(nu);
   else
    G.Tvu = mdsdata('STATIC("T_V"//$1//"_"//$2)[*,$3]',selu(2:end),selu,dimu);
   end
   G.Ru    = mdsdata('STATIC("RES_"//$1)[$2]',selu,dimu);
   G.Mfu   = mdsdata('STATIC("MUT_"//$1//"_"//$2)[$3,$4]' , self,selu,G.dimf,dimu);
   G.Ru    = mdsdata('STATIC("RES_"//$1)[$2]',selu,dimu);
   G.Mfu   = mdsdata('STATIC("MUT_"//$1//"_"//$2)[$3,$4]' , self,selu,G.dimf,dimu);
   G.Bmu   = mdsdata('STATIC("BT_M_"//$1) [$2,$3]' ,selu,G.dimm,dimu);
   
   % Retrieve matrix mapping vessel to segment currents to output Is
   % Is = Tivs*Iv; See JMM toolbox documentation
   G.Tivs = mdsdata('static("TINV_S_V")');
 end

 %% Optional extra output arguments only used by some codes
 for ii=1:numel(varargin)
   switch varargin{ii}
     case 'Maa'
       G.Maa   = mdsdata('STATIC("MUT_A_A"   )[ $1,$2]',G.dima,G.dima);
     case 'Ra'
       G.Ra    = mdsdata('STATIC("RES_A")[ $1]', G.dima);
     case 'Muu'
       G.Muu   = mdsdata('STATIC("MUT_"//$1//"_"//$2//"")[$3,$4]', selu,selu,dimu,dimu);
     case 'Mas'
       G.Mas   = mdsdata('STATIC("MUT_"//$1//"_S"       )[$1,$2]', G.dima,G.dims);
     case 'Mau'
       G.Mau   = mdsdata('STATIC("MUT_A_"//$1//""       )[$2,$3]', selu,G.dima,dimu);
     case 'Mxs'
        if ~isempty(selx) , G.Mxs   = mdsdata('STATIC("MUT_"//$1//"_S"       )[* ,$2]', selx,G.dims);   end
     case 'Mxu'
       if ~isempty(selx) && ~isempty(dimu), G.Mxu   = mdsdata('STATIC("MUT_"//$1//"_"//$2//"")[* ,$3]', selx,selu,dimu); end
     case 'Mxx'
       G.Mxx   = mdsdata('STATIC("MUT_"//$1//"_"//$1)       ',     selx       );
     case 'Brxa'
       if ~isempty(selx) && ~isempty(G.dima), G.Brxa  = mdsdata('STATIC("BR_"//$1//"_A"    )[ *,$2]',     selx,G.dima); end
     case 'Bzxa'
       if ~isempty(selx) && ~isempty(G.dima), G.Bzxa  = mdsdata('STATIC("BZ_"//$1//"_A"    )[ *,$2]',     selx,G.dima); end
     case 'Brxu'
        if ~isempty(selx) && ~isempty(dimu), G.Brxu  = mdsdata('STATIC("BR_"//$1//"_"//$2//"" )[* ,$3]', selx,selu,dimu); end
     case 'Bzxu'
        if ~isempty(selx) && ~isempty(dimu), G.Bzxu  = mdsdata('STATIC("BZ_"//$1//"_"//$2//"" )[* ,$3]', selx,selu,dimu); end
     case 'drBrxa'
        if ~isempty(selx) && ~isempty(G.dima), G.drBrxa = mdsdata('STATIC("DBRDR_"//$1//"_A"   )[* ,$2]', selx,G.dima); end
     case 'drBzxa'
        if ~isempty(selx) && ~isempty(G.dima), G.drBzxa = mdsdata('STATIC("DBZDR_"//$1//"_A"   )[* ,$2]', selx,G.dima); end
     case 'dzBrxa'
        if ~isempty(selx) && ~isempty(G.dima), G.dzBrxa = mdsdata('STATIC("DBRDZ_"//$1//"_A"   )[* ,$2]', selx,G.dima); end
     case 'dzBzxa'
        if ~isempty(selx) && ~isempty(G.dima), G.dzBzxa = mdsdata('STATIC("DBZDZ_"//$1//"_A"   )[* ,$2]', selx,G.dima); end        
     case 'drBrxu'
        if ~isempty(selx) && ~isempty(dimu), G.drBrxu = mdsdata('STATIC("DBRDR_"//$1//"_"//$2//"")[* ,$3]', selx,selu,dimu); end    
     case 'drBzxu'
        if ~isempty(selx) && ~isempty(dimu), G.drBzxu = mdsdata('STATIC("DBZDR_"//$1//"_"//$2//"")[* ,$3]', selx,selu,dimu); end   
     case 'dzBrxu'
        if ~isempty(selx) && ~isempty(dimu), G.dzBrxu = mdsdata('STATIC("DBRDZ_"//$1//"_"//$2//"")[* ,$3]', selx,selu,dimu); end   
     case 'dzBzxu'
        if ~isempty(selx) && ~isempty(dimu), G.dzBzxu = mdsdata('STATIC("DBZDZ_"//$1//"_"//$2//"")[* ,$3]', selx,selu,dimu); end   
     case 'Bms'
       G.Bms  = mdsdata('STATIC("BT_M_S" )[$1,$2]', G.dimm,G.dims);
     case 'Rs'
       G.Rs   = mdsdata('STATIC("RES_S")[$1]', G.dims);
     otherwise
       error('%s: requested variable %s not implemented',mfilename,varargin{ii})
   end
 end

 G.dimu = dimu;
 
 %% Limits
 % Minimum/maximum limits
 G.Ialim = mdsdata('STATIC("INOM_A")[$1]',G.dima); % Current limits
 G.Valim = mdsdata('STATIC("UNOM_A")[$1]',G.dima); % Voltage limits
 G.Talim = inf(numel(G.dima),1);
 G.Talim(contains(G.dima,'G_001')) = 2000.^2 * 0.23; % Thermal limit: 2000A during 0.23s (from Ugo Siravo discussion)
 
 iTOR=find(contains(G.dima,'TOR'));
 if ~isempty(iTOR)
   if G.Valim(iTOR)==0,   G.Valim(iTOR) = 1400; % manually set limit
   else, warning('G.Valim(iTOR) is no longer zero, someone fixed the static tree?'); % message for posterity
   end
 end
 
 if P.shot >= 100000 && P.shot < 200000 || P.shot < 25094
   % Fictitious shots only contain the PCS tree, not the BASE tree
   % Values for shots prior to 25094 are in \BASE::COIL_PROT:COEFF
   % and THRESHOLD. They contain 10 equations, scaling differs, some values
   % differ and matrix is transposed ...
   base_shot = -1;
 else
   base_shot = P.shot;
 end
 meqmdsopen(base_shot,'BASE',[],P.mdsserver);

 % Coil protection equation limits - see 'TCV Coil Protection' SPCwiki page
 [ceqdim,s] = mdsvalue('dim_of(\BASE::COIL_PROT:CEQ,1)');
 assert(isodd(s),'could not load coil protection equations for this shot %d',P.shot)
 % Identify coils without equations and equations with incomplete coil set
 [mask,iaeq] = ismember(deblank(ceqdim),G.dima);
 Ceq = mdsdata('\BASE::COIL_PROT:CEQ'); ieqok = find(~any(Ceq(:,~mask),2));
 G.Ceq(:,iaeq(mask))  = Ceq(ieqok,mask);G.Ceq(:,end+1:numel(G.dima)) = 0; % Pad remaining coils with 0s.
 % eliminate equations involving coils that are not present
 G.Seq  = mdsdata('\BASE::COIL_PROT:SEQ[$1]',ieqok);
 G.Ilow = mdsdata('\BASE::COIL_PROT:ILOW[$1]',ieqok);
 G.Ihig = mdsdata('\BASE::COIL_PROT:IHIG[$1]',ieqok);
 G.Ihig(isnan(G.Ihig)) = Inf;
end
