%LIUXTCV  TCV LIUQE diagnostic measurements
% LX = LIUXTCV(SHOT,T,L) returns a structure with measurements at times T
% for shots SHOT, taken from TCV_SHOT tree, or using LIUXRTC if L.P.XRTC. See
% also LIUX,LIUXRTC.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function X = liuxtcv(shot,t,L)
 if isfield(L.P,'xrtc') && L.P.xrtc
  X = liuxrtc(shot,t,L);
 else
  X.tokamak = 'TCV';
  % Default time selection
  if isempty(t)
   assert(isscalar(shot),'If multiple shots are provided, t cannot be empty');
   t = liuqe_times(shot);
   % Fallback if t is still empty
   if isempty(t)
    X.t = [];
    X.shot = [];
    return;
   end
  end
  if numel(shot) == 1, shot = repmat(shot,1,numel(t   )); end
  if numel(t   ) == 1, t    = repmat(t   ,1,numel(shot)); end
  nt = numel(t);
  shot = shot(:)';
  sl = fliplr(unique(shot));
  for s = sl
   ks = find(shot == s);
   % Open shot
   meqmdsopen(s,'TCV_SHOT',[],L.P.mdsserver)
   ts = t(ks);
   X.shot(ks) = s;
   % For TCV, all magnetics share the same time base.
   % Get time of measurements from first signal, then check consistency
   if L.P.reff == -1
    [X.Ff(1:L.G.nf,ks),txks] = mdsx('\MAGNETICS::RFLUX',ts,L.G.dimf);
    [X.Ff(1       ,ks)     ] = mdsx('\MAGNETICS::FLUX' ,ts,'001'   ,txks);
   else
    [X.Ff(1:L.G.nf,ks),txks] = mdsx('TCV_IDEALLOOP("FLUX" )',ts,L.G.dimf);
   end
   X.Bm(1:L.G.nm,ks)  = mdsx('\MAGNETICS::BPOL_003'  ,ts,L.G.dimm,txks);
   X.Ia(1:L.G.na,ks)  = mdsx('\MAGNETICS::IPOL'      ,ts,L.G.dima,txks);
   X.Uf(1:L.G.nf,ks)  = mdsx('TCV_IDEALLOOP("VLOOP")',ts,L.G.dimf,txks);
   X.rBt(1      ,ks)  = mdsx('\MAGNETICS::RBPHI'     ,ts,{}      ,txks);
   if L.P.idml == 0
    X.Ft(1,ks) = zeros(1,numel(ts));
   else
    X.Ft(1,ks) = mdsx('\RESULTS::DMLCOR',ts,{},txks) + ...
     mdsx(sprintf('\\RESULTS::DMLBO%03d',find(strcmp('FIX-PLA-LIN',deblank(mdsdata('GETNCI("\\RESULTS::DML:BO0%%.PARAM:TYPE","RECORD")'))))),...
     ts,{},txks);
   end
   % optional other measurements
   if L.P.kalman
     for ia = 1:L.G.na
       mydima = L.G.dima{ia};
       if startsWith(mydima,{'OH','E','F'})
         UDC1 = mdsx(sprintf('\\POWER::RED_U_DC1:%s',mydima),ts,{},txks);
         UDC2 = mdsx(sprintf('\\POWER::RED_U_DC2:%s',mydima),ts,{},txks);
         X.Va(ia,:) = UDC1+UDC2;
       elseif startsWith(mydima,'G')
         X.Va(ia,:) = mdsx('\POWER::FPS:UFOS',ts,{});
       end
     end
   end
   X.t(1,ks) = txks; % overwrite time with true magnetics time in time index
   if numel(sl) > 10, fprintf('LIUQE:LIUXTCV:Filled #%d (%d time slices)\n',s,numel(ts)), end
  end
  
  if ~L.P.ivesm
   X.Iu = zeros(L.G.nu,nt); % might replace by empty later?
  elseif isequal(L.P.selu(1),'s')
   X.Iu = X.Uf./(-L.G.Ru(:)); % Segment currents
  else
   error('Can''t use selu=%s with ivesm=1. Set selu=''s''/''sx'' or ivesm=0.',L.P.selu);
  end
  
  % For TCV-LIUQE Ip observer we always need the segment vessel currents
  % While we may have loaded another description in selu.
  assert(~isempty(L.G.Ipm),'must use Ip estimator for TCV');
  Is   = X.Uf./(-L.G.Rs(:)); 
  X.Ip =  L.G.Ipm*X.Bm + L.Ipa*X.Ia + L.Ips*Is;
 end 
end

function [x,tx] = mdsx(node,t,dim,txin)
 % Cache expression stored in node to avoid reading it again when populating T
 if isempty(dim), [X,stat] = mdsdata(['__MEQX=',node]);
 else,            [X,stat] = mdsdata(['__MEQX=',node '[*,$1]'],dim); end
 assert(isodd(stat) && ~isempty(X),'LIUQE:MDS:Cannot read measurement %s',node)
 [T,stat] = mdsdata('DIM_OF(__MEQX,0)');
 assert(~any(isnan(X(:))),'Invalid data for %s',node)
 assert(numel(T)==size(X,1),'LIUQE:MDS:Size mismatch between %s and its time dimension',node)
 
 assert(isodd(stat),'error loading DIM_OF(%s)',node);
 kt = min(max(round((t-T(1))/mean(diff(T))),0)+1,numel(T));
 x = NaN(size(X,2),numel(t));
 x(:,kt~=0) = X(kt(kt~=0),:)';
 % Time of measurement
 if nargout > 1 || nargin > 3
  tx = t;
  tx(kt~=0) = T(kt(kt~=0));
  % Check it matches previous signals
  if nargin > 3
   % DML results tree used to compute the time dimension slope, so could be slightly different
   % Magnetics are currently aquired at 10 kHz: dt=1e-4 >> 1e-6
   nomismatch =  (max(abs(txin-tx))<1e-6);
   if ~nomismatch && numel(t)>1 % check mean diagnostic sampling time w.r.t. what is requested
    Ts = mean(diff(T)); ts = min(diff(t));
    % Even if t and T are equal, the following condition can be violated
    % due to cumulative roundoff errors. But in this case, no time mismatch
    % should be found. So we only use this check as a possible hint for
    % reason behing the time mismatch.
    assert(round(Ts*1e6)<=round(ts*1e6),'liuxtcv:timeMismatch',...
      'Time array mismatch in liuxtcv for node %s. Sampling time is %dus, while asking for up to %dus',...
      node,round(Ts*1e6),round(ts*1e6));
   end
   assert(nomismatch,'liuxtcv:timeMismatch','Time array mismatch in liuxtcv for node %s',node);
  end
 end
end
