%LIUXCREATE CREATE(ITER) LIUQE diagnostic measurements
% LX = LIUXCREATE(EQUIL,T) returns a structure with measurements at time T
% for EQUIL(2:end), taken from CREATE equilibrium database. The first entry
% in EQUIL contains the sensors and coils geometric information, used by
% LIUGCREATE, LIUPCREATE when called in LIU.
% See also LIUX.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function LX = liuxcreate(equil,t,L)
% Locate equilibrium file (first entry in equil is the geometry file, see
% liupcreate, liugcreate).
createx = equil(2:end);
[file,equil,t] = createf(createx,t);
nt = numel(t);
if length(file)>1
  for k = nt:-1:1
    D = load(file{k});
    buildLX(k);
    LX.shot(:,k) = equil(k);
    LX.t   (:,k) = t(k);
  end
else
  D = load(file{1});
  k = 1 : nt;
  LX.shot = equil*ones(1,length(k));
  LX.t    = t;
  buildLX(k,t);
end


if not(isfield(D.coils_data, 'connection'))
  % Remove CS1L and VS3L (default for ITER)
  k = [find(strcmp(D.coils_data.names,'coil CS1L')) ;
    find(strcmp(D.coils_data.names,'coil VS3L'))];
else
  k = [];
  for i = 1 : size(D.coils_data.connection,1)
    k = [k; find(strcmp(D.coils_data.names,['coil ' D.coils_data.connection{i,1}]))];
  end
end
LX.Ia(k,:) = [];


  function buildLX(k,t)
    
    if nargin>1 && isfield(D,'time')
      if numel(D.time) == numel(t)
        tidx = 1:numel(t);
      else
        tidx = find(D.time>=t(1) & D.time<=t(end));
      end
    else
      tidx = 1:nt;
    end
    
    if ~iscell(D.saddle(1).name) % for retrocompatibility
      for i = 1 : numel(D.saddle), D.saddle(i).name = {D.saddle(i).name}; end
    end
    ksel = find(~cellfun(@isempty,regexp([D.saddle.name],L.P.selfm)));
    for kk = 1:length(ksel)
      saddle_values(kk,:) = D.saddle(ksel(kk)).values(tidx);
    end
    %   LX.Ff  (:,k) = reshape(saddle_values,length(k),[])';
    LX.Ff(:,k) = saddle_values;
    ksel = find(~cellfun(@isempty,regexp(D.field_sensors.names,L.P.selfm)));
    LX.Bm  (:,k) = reshape( D.field_sensors.values(ksel,tidx),[],length(k));
    LX.Ia  (:,k) = reshape( D.coils_data.current(:,tidx)     ,[],length(k));
    LX.rBt (:,k) =          D.Bvacradius(tidx)                ;
    if isfield(D.diamagneticflux,'values')
      LX.Ft  (:,k) =        D.diamagneticflux.values(1,tidx)  ;
    else
      LX.Ft  (:,k) =    abs(D.diamagneticflux(tidx))          ;   % TODO: temporary fix for ITER case 3, cases 1,5,7,8 have all positive fluxes
    end
    LX.Is  (:,k) =                                  zeros(0,length(k));
    LX.Ip  (:,k) =          D.plasmacurrent(tidx)             ;
    try
      LX.li  (:,k) =          D.globalParameters.li3ITER.value(tidx);
      LX.bp  (:,k) =          D.globalParameters.betap(tidx)      ;
      LX.rA  (:,k) =          D.globalParameters.magneticAxis(1,tidx);
      LX.zA  (:,k) =          D.globalParameters.magneticAxis(2,tidx);
      LX.PpQ (:,k) = reshape( D.fluxFunctionProfiles.staticPprime(:,tidx)          ,[],length(k));
      LX.TTpQ(:,k) = reshape( D.fluxFunctionProfiles.ffprime(:,tidx)               ,[],length(k));
      LX.FNQ (:,k) = reshape( D.fluxFunctionProfiles.normalizedPoloidalFlux(:,tidx),[],length(k));
      %    figure('visible','off')
      %    h = pdecont(D.meshandequil.p,D.meshandequil.t,D.meshandequil.psi,[D.meshandequil.psb D.meshandequil.psb]);
      %    LX.rc(1:length(h(1).XData),k) = h(1).XData;
      %    LX.zc(1:length(h(1).YData),k) = h(1).YData;
      %    close(gcf)
    end
    
    % no vessel current mesurements
    LX.Iu = zeros(L.G.nu,numel(LX.t));
    
  end
end


