%LIUXRTC  TCV LIUQE diagnostic measurements from RTC tree
% LX=liuxrtc(shot,t,L) loads LIUQE input data from RTC tree, 
% for specific shots in which data was stored in mems.
% If t is empty, returns the data on the same time basis as stored in the tree
%
%  * For SCD Mark 2 system: (2019-..):
%     Parses all MEM names for occurence of 'LX'
%  * For SCD Mark 1 system: (2009-2019):
%     Assumes shot was run with node03 running model that stores LIUQE inputs
%     on mems. Searches for 'liuqe_inp' or 'LX' in mems name.
%
% EXAMPLES: 
% Mark 2 system:
%  Get data from mems of RT02 thread 03.
%  LX = liuxrtc(65488,t,L);
%
% Mark 1 system (shot<65300):
%  Get data from mems of RT03.
%   LX = liuxrtc(57587,t,L);
%
% See also TS2MEQ
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function LX = liuxrtc(shot,t,L)
narginchk(1,3)

if shot < 65300
  %% get data from mems - old SCD
  LX = liuxscdmems_pre2019(shot,L);
else
  LX = liuxscdmems(shot,L);
end

% interpolate if necessary
if ~isempty(t)
  LX = meqinterp(LX,t,'nearest');
end

LX.tokamak = 'TCV';

end

function LX = liuxscdmems(shot,L)
  fprintf('get LX data from RTC tree stored MEMs for shot %d\n',shot);
  assert(~isempty(which('get_scd_mems')),'liuxrtc:noscd','did not find get_scd_mems.m - need to add SCD tools to the path')
  LXts = getfield(get_scd_mems(shot,'LX','verbosity',0),'LX'); % silent
  LX = ts2meq(L,LXts,[]);
  LX.shot = shot;
  
  % split Ia,Iu if necessary
  if isfield(LX,'Ie')
    na = L.G.na; nu = L.G.nu;
    assert(size(LX.Ie,1)==na+nu,'incorrect size for LX.Ie')
    LX.Ia = LX.Ie(1:na    ,:);
    LX.Iu = LX.Ie(na+1:end,:);
  elseif isfield(LX,'Is')
    LX.Iu = LX.Is;
  end
end

function LX = liuxscdmems_pre2019(shot,L)

fprintf('get LX data for shot %d from RTC tree stored MEMs for CRPPRT03\n',shot);

meqmdsopen(shot,'rtc',[],'SCD');

% get time steps
t = mdsvalue('dim_of(\top.crpprt03.thread1.mems.mem_001,0)');
nt = numel(t);

% Identify LIUQE inputs in mems
[names,status] = mdsvalue('getnci("\\top.crpprt03.thread1.mems:mem_%%%:name","record")');
assert(isodd(status),'could not get mems names');
% Assume we get all of them in order (safe) so first name corresponds to mem_001 etc.
patternsLiuqeInput = {'liuqe_inp','LX'};
flag = contains(names,patternsLiuqeInput);
ifirstLX = find(flag,1,'first');
assert(~isempty(ifirstLX),'could not find liuqe inputs in mems by parsing name strings - shot was probably not run with correct code on node03');
nInputs = find(flag,1,'last') - ifirstLX + 1;

% there are two cases: if the 'old' storage way is used, then only 18 PF
% coils are stored, and there are 134 inputs.
if nInputs == 134
  gCoilAvailable = false;
elseif nInputs == 135
  gCoilAvailable = true;
else
  error('unknown or strange number of LIUQE inputs detected in mem');
end

% Ff
istart = ifirstLX-1;
n = L.G.nf; % number of Flux loops
Ff = zeros(n,nt);
for ii=1:n
   imem = istart + ii;
   Ff(ii, :) = mdsvalue(sprintf('\\top.crpprt03.thread1.mems.mem_%03d',imem));
end

% Bm
istart = istart + 38;
n = L.G.nm;
Bm = zeros(n,nt);
for ii=1:n
   imem = istart + ii;
   Bm(ii, :) = mdsvalue(sprintf('\\top.crpprt03.thread1.mems.mem_%03d',imem));
end

% Ia
istart = istart + 38;
if gCoilAvailable
 n = 19;
else
 n = 18;
end

Ia = zeros(n,nt);
for ii=1:n
   imem = istart + ii;
   Ia(ii, :) = mdsvalue(sprintf('\\top.crpprt03.thread1.mems.mem_%03d',imem));
end

if ~gCoilAvailable
  warning('G coil data was not avaialble in mems for this shot - substituting with zeros');
  % add vector with zeros for Ia
  Ia = [Ia(1:16,:);zeros(1,size(Ia,2));Ia([17,18],:)];
end

% Is
istart = istart + n;
n = 38;
Is = zeros(n,nt);
for ii=1:n
   imem = istart + ii;
   Is(ii, :) = mdsvalue(sprintf('\\top.crpprt03.thread1.mems.mem_%03d',imem));
end

% Ft
istart = istart + n;
n = 1;
Ft = zeros(n,nt);
for ii=1:n
   imem = istart + ii;
   Ft(ii, :) = mdsvalue(sprintf('\\top.crpprt03.thread1.mems.mem_%03d',imem));
end

% rBt
istart = istart + n;
n = 1;
rBt = zeros(n,nt);
for ii=1:n
   imem = istart + ii;
   rBt(ii, :) = mdsvalue(sprintf('\\top.crpprt03.thread1.mems.mem_%03d',imem));
end

%% Update to most recent LX style
Iu = Is;
assert(strcmp(L.P.selu(1),'s'),'must use selu=''s'' for loading LX data from rtc mems')

Ip = L.G.Ipm*Bm + L.Ipa*Ia + L.Ips*Is;
Uf = Iu.*(-L.G.Ru(:)); % Segment currents

%% save all in LX structure
LX = struct('tokamak','TCV','shot',shot,'t',t);
LX = meqlarg(LX,Ff,Bm,Ia,Iu,Uf,Ft,rBt,Ft,Ip);

mdsclose;
mdsdisconnect;

end
