gdat_imas.m

function [gdat_data,gdat_params,error_status,varargout] = gdat_imas(shot,data_request,varargin)
%
% function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
%
% Aim: get data from a given machine using full path or keywords.
%      data_request are and should be case independent (transformed in lower case in the function and outputs)
%
% If no inputs are provided, return the list of available pre-defined data_request in gdat_data and default parameters gdat_params
%
% IMAS means assume access data from IDS database ala IMAS with tokamakname and dataversion and run_number using imas_open_env and ids_get calls
% Should be the same on ITER hpc and gateway, to be checked
%
% At first only "ids" request implemented, but keywords like ip etc should relate to specific IDS as for any machine
%
% Inputs:
%
%    no inputs: return default parameters in a structure form in gdat_params
%    shot: shot number (if empty uses 'idx_imas_open_env' in params from previous call or input param)
%    data_request: keyword (like 'ip') or trace name or structure containing all parameters but shot number
%    varargin{i},varargin{i+1},i=1:nargin-2: additional optional parameters given in pairs: param_name, param_value
%                                            The optional parameters list might depend on the data_request
%              examples of optional parameters:
%                               'plot',1 (plot is set by default to 0)
%                               'machine','IMAS' (the default machine is the local machine)
%
%
% Outputs:
%
% gdat_data: structure containing the data, the time trace if known and other useful information
% gdat_data.t : time trace
% gdat_data.data: requested data values
% gdat_data.dim : values of the various coordinates related to the dimensions of .data(:,:,...)
%                     note that one of the dim is the time, replicated in .t for clarity
% gdat_data.dimunits : units of the various dimensions, 'dimensionless' if dimensionless
% gdat_data.error_bar : if provided
% gdat_data.gdat_call : list of parameters provided in the gdat call (so can be reproduced)
% gdat_data.shot: shot number
% gdat_data.machine: machine providing the data
% gdat_data.gdat_request: keyword for gdat if relevant
% gdat_data.data_fullpath: full path to the data node if known and relevant, or relevant expression called if relevant
% gdat_data.gdat_params: copy gdat_params for completeness (gdat_params contains a .help structure detailing each parameter)
% gdat_data.xxx: any other relevant information
%
% gdat_params contains the options relevant for the called data_request. It also contains a help structure for each option
% eg.: param1 in gdat_params.param1 and help in gdat_params.help.param1
%
% Note for liuqe parameter: At this time we have liuqe_fortran (which used to be the default) and liuqe_matlab (which is new and now becomes/is the default)
% We still have liuqe1, 2 and 3 nodes for either of them, so you can choose:
% 'liuqe',1 (2 or 3): to get the default liuqe 1, 2 or 3 (which is now the matlab version)
% 'liuqe',11 (12 or 13): to get liuqe_fortran nodes 1, 2 or 3
% 'liuqe',21 (22 or 23): to get liuqe_matlab nodes 1, 2 or 3 (may be needed if we set the default to liuqe_fortran for old shots)
% 'liuqe',-1  : to get FBTE result
%
% Examples:
% (should add working examples for various machines (provides working shot numbers for each machine...))
%
%    to get the list of ids: aa=gdat([],'ids','machine','imas'); aa.gdat_params.sources_available
%    to get an empty ids:    aa=gdat([],'ids','machine','imas','source','core_profiles'); aa.core_profiles
%
%    [a1,a2]=gdat;
%    a2.data_request = 'ids'; a2.source='equilibrium'; a2.machine = 'imas';
%    a3=gdat(48836,a2);  % gives input parameters as a structure, allows to call the same for many shots
%    a4=gdat('opt1',123,'opt2',[1 2 3],'shot',48832,'data_request','Ip','opt3','aAdB'); % all in pairs
%    a5=gdat(48836,'ids','source','core_profiles'); % standard call for a given IDS core_profiles, use {'',''} for a list of IDSs
%    a6=gdat(48836,'ip'); % standard call for a specific keyword
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% Comments for local developer:
% gdat is just a "header routine" calling the gdat for the specific machine gdat_`machine`.m which can be called
% directly, thus which should be able to treat the same type of input arguments
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prepare some variables, etc

varargout{1}=cell(1,1);
error_status=1;

% construct main default parameters structure
gdat_params.data_request = '';
default_machine = 'imas';

gdat_params.machine=default_machine;
gdat_params.doplot = 0;
gdat_params.nverbose = 1;
gdat_params.run = []; % should be defined by user to make sure
gdat_params.occurence = 0; % default/usual value

% special case need to treat here if params given as structure in data_request:
if isstruct(data_request)
  gdat_params = data_request;
end

% construct list of keywords from global set of keywords and specific IMAS set
% get data_request names from centralized table
%%% data_request_names = get_data_request_names; % do not use xlsx anymore but scan respective machine_requests_mapping.m files
data_request_names = get_data_request_names_from_gdat_xxx(default_machine);
% add IMAS specific to all:
if ~isempty(data_request_names.imas)
  imas_names = fieldnames(data_request_names.imas);
  for i=1:length(imas_names)
    data_request_names.all.(imas_names{i}) = data_request_names.imas.(imas_names{i});
  end
end
data_request_names_all = sort(fieldnames(data_request_names.all));

% construct default output structure
gdat_data.data = [];
gdat_data.units = [];
gdat_data.dim = [];
gdat_data.dimunits = [];
gdat_data.t = [];
gdat_data.x = [];
gdat_data.shot = [];
gdat_data.gdat_request = [];
gdat_data.gdat_params = gdat_params;
gdat_data.data_fullpath = [];


% Treat inputs:
ivarargin_first_char = 3;
data_request_eff = '';
if nargin>=2
  if ischar(data_request);
    data_request_eff = data_request;
    data_request = data_request; % should not lower until see if keyword
  elseif isstruct(data_request)
    data_request_eff = data_request.data_request;
  end
end
gdat_data.gdat_request = data_request_names_all; % so if return early gives list of possible request names
gdat_data.gdat_params.help = imas_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_params = gdat_data.gdat_params;

% no inputs
if nargin==0
  % return defaults and list of keywords
  return
end

do_mdsopen_mdsclose = 1;
% treat 1st arg
if nargin>=1
  if isempty(shot)
    % if shot already open, will use idx from gdat_params, otherwise returns
    do_mdsopen_mdsclose = 0;
    if ~isfield(gdat_data.gdat_params,'idx_imas_open_env')
      if ~strcmp(data_request_eff,'ids'); return; end % empty shot return empty ids so valid command
    end
  elseif isnumeric(shot)
    gdat_data.shot = shot;
  elseif ischar(shot)
    ivarargin_first_char = 1;
  else
    if gdat_params.nverbose>=1; warning('type of 1st argument unexpected, should be numeric or char'); end
    error_status=2;
    return
  end
  if nargin==1
    % Only shot number given. If there is a default data_request set it and continue, otherwise return
    return
  end
end
% 2nd input argument if not part of pairs
if nargin>=2 && ivarargin_first_char~=1
  if isempty(data_request)
    return
  end
  % 2nd arg can be a structure with all options except shot_number, or a string for the pathname or keyword, or the start of pairs string/value for the parameters
  if isstruct(data_request)
    if ~isfield(data_request,'data_request')
      if gdat_params.nverbose>=1; warning('expects field data_request in input parameters structure'); end
      error_status=3;
      return
    end
    data_request_eff = data_request.data_request;
    % some request are case sensitive because mds names can be like \magnetics::ipol[*,"OH_001"]
    data_request.data_request = data_request.data_request;
    gdat_params = data_request;
  else
    % since data_request is char check from nb of args if it is data_request or a start of pairs
    if mod(nargin-1,2)==0
      ivarargin_first_char = 2;
    else
      ivarargin_first_char = 3;
      data_request_eff = data_request;
    end
  end
end

if ~isstruct(data_request)
  gdat_params.data_request = data_request_eff;
end

% if start pairs from shot or data_request, shift varargin
if ivarargin_first_char==1
  varargin_eff{1} = shot;
  varargin_eff{2} = data_request;
  varargin_eff(3:nargin) = varargin(:);
elseif ivarargin_first_char==2
  varargin_eff{1} = data_request;
  varargin_eff(2:nargin-1) = varargin(:);
else
  varargin_eff(1:nargin-2) = varargin(:);
end

% extract parameters from pairs of varargin:
if (nargin>=ivarargin_first_char)
  if mod(nargin-ivarargin_first_char+1,2)==0
    for i=1:2:nargin-ivarargin_first_char+1
      if ischar(varargin_eff{i})
	gdat_params.(lower(varargin_eff{i})) = varargin_eff{i+1}; % cannot enforce lower(params value) if filename or else are case sensitive
      else
        if gdat_params.nverbose>=1; warning(['input argument nb: ' num2str(i) ' is incorrect, expects a character string']); end
        error_status=401;
        return
      end
    end
  else
    if gdat_params.nverbose>=1; warning('number of input arguments incorrect, cannot make pairs of parameters'); end
    error_status=402;
    return
  end
end

% if it is a request_keyword copy it:
if ischar(data_request_eff) || length(data_request_eff)==1
  ij=strmatch(lower(data_request_eff),data_request_names_all,'exact');
else
  ij=[];
end
if ~isempty(ij);
  gdat_data.gdat_request = data_request_names_all{ij};
  gdat_params.data_request = gdat_data.gdat_request;
  if isfield(data_request_names.all.(data_request_names_all{ij}),'description') && ~isempty(data_request_names.all.(data_request_names_all{ij}).description)
    % copy description of keyword
    gdat_data.request_description = data_request_names.all.(data_request_names_all{ij}).description;
  end
end

% special treatment if shot and data_request given within pairs
if isfield(gdat_params,'shot')
  shot = gdat_params.shot; % should use only gdat_params.shot but change shot to make sure
  gdat_data.shot = gdat_params.shot;
  gdat_params=rmfield(gdat_params,'shot');
end

if ~isfield(gdat_params,'data_request') || isempty(gdat_params.data_request)
  % warning('input for ''data_request'' missing from input arguments') % might be correct, asking for list of requests
  error_status=5;
  return
end
gdat_data.gdat_params = gdat_params; % This means that before here it is gdat_params which should be updated
clear gdat_params; % to make sure one uses gdat_data.gdat_params now

% re-assign main variables to make sure use the one in gdat_data structure
shot = gdat_data.shot;
data_request_eff = gdat_data.gdat_params.data_request;
if ~isfield(gdat_data.gdat_params,'database_user')
  gdat_data.gdat_params.database_user = getenv('USER');
  if gdat_data.gdat_params.nverbose >= 3; disp('set database_user to username as default since not provided'); end
end
if ~isfield(gdat_data.gdat_params,'tokamak_name')
  gdat_data.gdat_params.tokamak_name = gdat_data.gdat_params.machine;
  if gdat_data.gdat_params.nverbose >= 3; disp('set tokamak_name to machine as default since not provided'); end
end
if ~isfield(gdat_data.gdat_params,'data_major_version')
  gdat_data.gdat_params.data_major_version = '3';
  if gdat_data.gdat_params.nverbose >= 3; disp('set data_major_version to ''3'' as default since not provided'); end
end
error_status = 6; % at least reached this level
gdat_params = gdat_data.gdat_params;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Specifications on how to get the data provided in imas_requests_mapping
mapping_for_imas = imas_requests_mapping(data_request_eff,shot);
gdat_data.label = mapping_for_imas.label;
ishot=NaN;
if do_mdsopen_mdsclose
  if ~isfield(gdat_data.gdat_params,'run') || isempty(gdat_data.gdat_params.run)
    warning('run number not provided, cannot continue')
    return
  end
  try
    gdat_data.gdat_params.idx_imas_open_env = imas_open_env('ids', shot, gdat_data.gdat_params.run,gdat_data.gdat_params.database_user, ...
	  gdat_data.gdat_params.tokamak_name,gdat_data.gdat_params.data_major_version);
    disp('')
  catch ME
    warning('could not imas_open_env with following gdat_data.gdat_params:');
    gdat_data.gdat_params
    rethrow(ME)
  end
  if gdat_data.gdat_params.idx_imas_open_env < 0
    if gdat_data.gdat_params.nverbose>=1
      warning(['idx < 0, shot/run not opened properly with gdat_data.gdat_params:']);
      gdat_data.gdat_params
    end
    return
  end
else
  if isfield(gdat_data.gdat_params,'idx_imas_open_env')
    if gdat_data.gdat_params.idx_imas_open_env < 0
      if gdat_data.gdat_params.nverbose>=1
	warning(['idx < 0, shot/run not opened properly with gdat_data.gdat_params:']);
	gdat_data.gdat_params
      end
      return
    end
  else
    if ~isempty(shot)
      if gdat_data.gdat_params.nverbose>=1
	warning(['no idx in gdat_params, cannot open shot with gdat_data.gdat_params:']);
	gdat_data.gdat_params
      end
      return
    end
  end
end
gdat_params = gdat_data.gdat_params;

% fill again at end to have full case, but here to have present status in case of early return
gdat_data.gdat_params.help = imas_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_data.mapping_for.imas = mapping_for_imas;
gdat_params = gdat_data.gdat_params;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(mapping_for_imas.method,'expression')
  % method="expression"
  % assume expression contains an expression to evaluate and which creates a local structure into variable gdat_tmp
  % we copy the structure, to make sure default nodes are defined and to avoid if return is an closed object like tdi
  % eval_expr = [mapping_for_imas.expression ';'];
  eval([mapping_for_imas.expression ';']);
  if isempty(gdat_tmp) || (~isstruct(gdat_tmp) & ~isobject(gdat_tmp))
    if (gdat_params.nverbose>=1); warning(['expression does not create a gdat_tmp structure: ' mapping_for_imas.expression]); end
    error_status=801;
    return
  end
  tmp_fieldnames = setdiff(fieldnames(gdat_tmp),{'gdat_request','label'}); % could/should also remove label in any case
  if sum(strcmp(tmp_fieldnames,'data'))==0 % note: cannot do isfield since gdat_tmp might be an object
    if (gdat_params.nverbose>=1); warning(['expression does not return a child name ''data'' for ' data_request_eff]); end
  end
  for i=1:length(tmp_fieldnames)
    gdat_data.(tmp_fieldnames{i}) = gdat_tmp.(tmp_fieldnames{i});
  end
  % add .t and .x in case only dim is provided
  % do not allow shifting of timedim since should be treated in the relevant expression
  ijdim=find(strcmp(tmp_fieldnames,'dim')==1);
  if ~isempty(ijdim)
    nbdims = length(gdat_data.dim);
    if mapping_for_imas.timedim==-1;
      mapping_for_imas.timedim = nbdims;
      if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_imas.timedim = nbdims-1; end
    end
    dim_nontim = setdiff([1:nbdims],mapping_for_imas.timedim);
    ijt=find(strcmp(tmp_fieldnames,'t')==1);
    if isempty(ijt)
      gdat_data.t = gdat_data.dim{mapping_for_imas.timedim};
    end
    ijx=find(strcmp(tmp_fieldnames,'x')==1);
    if isempty(ijx)
      if ~isempty(dim_nontim)
        % since most cases have at most 2d, copy as array if data is 2D and as cell if 3D or more
        if length(dim_nontim)==1
          gdat_data.x = gdat_data.dim{dim_nontim(1)};
        else
          gdat_data.x = gdat_data.dim(dim_nontim);
        end
      end
    end
    gdat_data.data_fullpath=mapping_for_imas.expression;
    if isfield(gdat_tmp,'help')
      gdat_data.help = gdat_tmp.help;
    else
      gdat_data.help = [];
    end
    if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
      [gdat_data] = gdat2time_out(gdat_data,1);
    end
  end
  % end of method "expression"

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_imas.method,'switchcase')
  switch data_request_eff % not lower(...) since data_request_eff should be lower case already at this stage
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % First the request names valid for "all" machines:
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'ids'}
    ids_empty_path = fullfile(fileparts(mfilename('fullpath')),'..','IMAS_IMAS','ids_empty');

    params_eff = gdat_data.gdat_params;
    if isfield(params_eff,'source') && ~isempty(params_eff.source)
      ids_top_name = params_eff.source;
    else
      ids_top_name = [];
      warning('gdat:EmptyIDSName','Need an ids name in ''source'' parameter\n check substructure gdat_params.sources_available for an ids list');
      addpath(ids_empty_path);
      assert(~~exist('ids_list','file'),'could not find ids_list.m in %s',ids_empty_path);
      gdat_data.gdat_params.sources_available = ids_list;
      rmpath(ids_empty_path);
      return
    end
    ids_gen_ok = ~~exist('ids_gen','file');

    if ids_gen_ok
      ids_empty = ids_gen(ids_top_name); % generate ids from gateway function ids_gen
    else
      % load empty ids structure from template file
      fname = sprintf('ids_empty_%s',ids_top_name);

      try
        assert(~~exist(ids_empty_path,'dir'),'folder %s not found',ids_empty_path);
        addpath(ids_empty_path);
        assert(~~exist(fname,'file'),'file %s does not exist in %s',fname,ids_empty_path);

        ids_empty = eval(fname);
        disp(['use structure in .mat file: ' ids_empty_path '/' fname]);
        rmpath(ids_empty_path);
      catch ME
        fprintf('Could not load empty template for %s\n',ids_top_name);
        fprintf('Available templates:\n');
        disp(dir(ids_empty_path));
        rmpath(ids_empty_path);
        rethrow(ME);
      end
    end
    if ~isfield(gdat_data.gdat_params,'error_bar') || isempty(gdat_data.gdat_params.error_bar)
      gdat_data.gdat_params.error_bar = 'delta';
    end
    if ~isfield(gdat_data.gdat_params,'cocos_in') || isempty(gdat_data.gdat_params.cocos_in)
      gdat_data.gdat_params.cocos_in = 11;
    end
    if ~isfield(gdat_data.gdat_params,'cocos_out') || isempty(gdat_data.gdat_params.cocos_out)
      gdat_data.gdat_params.cocos_out = 11;
    end
    if ~isfield(gdat_data.gdat_params,'ipsign_out') || isempty(gdat_data.gdat_params.ipsign_out)
      gdat_data.gdat_params.ipsign_out = 0;
    end
    if ~isfield(gdat_data.gdat_params,'b0sign_out') || isempty(gdat_data.gdat_params.b0sign_out)
      gdat_data.gdat_params.b0sign_out = 0;
    end
    if ~isfield(gdat_data.gdat_params,'ipsign_in') || isempty(gdat_data.gdat_params.ipsign_in)
      if gdat_data.gdat_params.ipsign_out~=0
        gdat_data.gdat_params.ipsign_in = mdsvalue('\pcs::data:iohfb');
      else
        gdat_data.gdat_params.ipsign_in = 0;
      end
    end
    if ~isfield(gdat_data.gdat_params,'b0sign_in') || isempty(gdat_data.gdat_params.b0sign_in)
      if gdat_data.gdat_params.b0sign_out~=0
        gdat_data.gdat_params.b0sign_in = mdsvalue('\pcs::data:if36fb');
      else
        gdat_data.gdat_params.b0sign_in = 0;
      end
    end
    try
      if isfield(gdat_data.gdat_params,'idx_imas_open_env') && gdat_data.gdat_params.idx_imas_open_env >= 0
	if gdat_data.gdat_params.occurence == 0
	  ids_top = ids_get(gdat_data.gdat_params.idx_imas_open_env,ids_top_name);
	else
	  ids_top = ids_get(gdat_data.gdat_params.idx_imas_open_env,[ids_top_name '/' num2str(gdat_data.gdat_params.occurence)]);
	end
	gdat_data.(ids_top_name) = ids_top;
      else
        gdat_data.(ids_top_name) = ids_empty;
      end
    catch ME_imas_ids_get
      disp(['there is a problem with: imas_get_ids_' ids_top_name ...
            ' , may be check if it exists in your path or test it by itself'])
      gdat_data.(ids_top_name) = ids_empty;
      gdat_data.([ids_top_name '_help']) = getReport(ME_imas_ids_get);
      rethrow(ME_imas_ids_get)
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   otherwise
    if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_imas']); end
    error_status=901;
    return
  end

else
  if (gdat_params.nverbose>=1); warning(['IMAS method=' mapping_for_imas.method ' not known yet, contact Olivier.Sauter@epfl.ch']); end
  error_status=602;
  return
end

gdat_data.gdat_params.help = imas_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_data.mapping_for.imas = mapping_for_imas;
gdat_params = gdat_data.gdat_params;
error_status=0;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [firthomratio] = get_fir_thom_rat_data(shot,maintracename,timebase,nverbose);
%
% since depends on shot number for using auto fit and thomson or thomson edge, use tracename and function here
%
% maintracename = 'thomson' or 'thomson_edge
%
% return fir_to_thomson ratio on time=timebase
%
% normally should use "main" thomson one built from auto fit if possible
% take edge one if need be
%
% default: maintracename = "thomson"
%
maintracename_eff = 'thomson';
if exist('maintracename') && ~isempty(maintracename)
  maintracename_eff = maintracename;
end

firthomratio = NaN;
if ~exist('timebase') || isempty(timebase)
  if (nverbose>=1)
    disp('need a timebase in get_fir_thom_rat_data')
  end
  return
end
firthomratio = NaN*ones(size(timebase));

if strcmp(maintracename_eff,'thomson')
  if shot>=23801
    tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
    if isempty(tracefirrat.data) || ischar(tracefirrat.data)
      if (nverbose>=1); disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat'); end
      tracefirrat=tdi('\results::thomson:fir_thom_rat');
    end
  else
    tracefirrat=tdi('\results::thomson:fir_thom_rat');
    if isempty(tracefirrat.data) || ischar(tracefirrat.data)
      if (nverbose>=1); disp('problem with \results::thomson:fir_thom_rat: empty'); end
    end
  end
elseif strcmp(maintracename_eff,'thomson_edge')
  if shot>=23801
    tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
    if isempty(tracefirrat.data) || ischar(tracefirrat.data)
      if (nverbose>=1); disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat'); end
      tracefirrat=tdi('\results::thomson:fir_thom_rat');
    end
  else
    tracefirrat=tdi('\results::thomson_edge:fir_thom_rat');
    if isempty(tracefirrat.data) || ischar(tracefirrat.data)
      if (nverbose>=1); disp('problem with \results::thomson_edge:fir_thom_rat: empty'); end
    end
  end
else
  if (nverbose>=1); disp('bad input in get_fir_thom_rat_data'); end
  return
end

if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) && any(isfinite(tracefirrat.data)) ...
      && ~isempty(tracefirrat.dim) && ~isempty(tracefirrat.dim{1})
  firthomratio = interp1(tracefirrat.dim{1},tracefirrat.data,timebase);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,doedge,nverbose);
%
try
  time=mdsdata('\results::thomson:times');
catch
  gdat_data.error_bar = [];
  if strcmp(data_request_eff(1:2),'ne')
    tracefirrat_data = [];
    gdat_data.firrat=tracefirrat_data;
    gdat_data.data_raw = gdat_data.data;
    gdat_data.error_bar_raw = gdat_data.error_bar;
  end

  if (nverbose>=1) && shot<100000
    warning('Problems with \results::thomson:times')
    warning(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
  end
  return
end
if isempty(time) || ischar(time)
  thomsontimes=time;
  if (nverbose>=1) && shot<100000
    warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
    disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
  end
  return
end
edge_str_ = '';
edge_str_dot = '';
if doedge
  edge_str_ = '_edge';
  edge_str_dot = '.edge';
end

nodenameeff = ['\results::thomson' edge_str_dot ':' data_request_eff(1:2)];
tracetdi=tdi(nodenameeff);
if isempty(tracetdi.data) || ischar(tracetdi.data) || isempty(tracetdi.dim)
  gdat_data.error_bar = [];
  gdat_data.firrat = [];
  gdat_data.data_raw = [];
  gdat_data.error_bar_raw = [];
  return
end

gdat_data.data=tracetdi.data'; % Thomson data as (t,z)
tracestd=tdi(['\results::thomson'  edge_str_dot ':' data_request_eff(1:2) ':error_bar']);
gdat_data.error_bar=tracestd.data';
gdat_data.data_fullpath=[nodenameeff '; error_bar'];
% add fir if ne requested
if strcmp(data_request_eff(1:2),'ne')
  tracefirrat_data = get_fir_thom_rat_data(shot,['thomson' edge_str_],time,nverbose);
  gdat_data.firrat=tracefirrat_data;
  gdat_data.data_raw = gdat_data.data;
  gdat_data.data = gdat_data.data_raw * diag(tracefirrat_data);
  gdat_data.error_bar_raw = gdat_data.error_bar;
  gdat_data.error_bar = gdat_data.error_bar_raw * diag(tracefirrat_data);
  gdat_data.data_fullpath=[gdat_data.data_fullpath '; fir_thom_rat ; _raw without firrat'];
  ij=find(isfinite(tracefirrat_data));
  if isempty(ij)
    gdat_data.data_fullpath=[gdat_data.data_fullpath ' WARNING use ne_raw since firrat has NaNs thus ne=NaNs; fir_thom_rat ; _raw without firrat'];
    disp('***********************************************************************')
    disp('WARNING: ne from Thomson has fir_thom_rat with Nans so only ne_raw can be used');
    disp('***********************************************************************')
  end
end
z=mdsdata(['\diagz::thomson_set_up' edge_str_dot ':vertical_pos']);
gdat_data.dim=[{z};{time}];
gdat_data.dimunits=[{'Z [m]'} ; {'time [s]'}];
gdat_data.x=z;
gdat_data.t=time;
% isfield does not work since tracetdi is not a 'struct' but a tdi object, thus use fieldnames
if any(strcmp(fieldnames(tracetdi),'units'))
  gdat_data.units=tracetdi.units;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, system_str_dot, psitbx_str, gdat_params )
% Input arguments:
%   gdat_data: current structure containing TS data
%   zshift: vertical shift of the equilibrium, can be a vector for time-dependent data
%   system_str_dot: '' or '.edge' indicating thomson syste,
%   psitbx_str: source for the equilibrium map, psitbx-style
%   gdat_params: parameters for the gdat call % TODO: Why is it different from gdat_data.gdat_params?

% Use psitbx
t_th = gdat_data.t;
th_z = mdsdata(['dim_of(\thomson' system_str_dot ':te,1)']); % TODO: Isn't this gdat_data.dim{1}?
th_r = 0.9*ones(size(th_z));
ntt = numel(t_th);
nch = numel(th_z); % number of chords

[t_psi,status] = mdsvalue('dim_of(imas_eq("i_p",$))',psitbx_str); % TODO: replace with time_psi when imas_eq supports it
if ~rem(status,2)
  warning('problem with getting equilibrium time base in gdat_imas')
  display(t_psi);
  return
end

zshifteff=zshift;
% set zshifteff time array same as t_psi
switch numel(zshifteff)
  case 1
    zshifteff=zshifteff * ones(size(t_th));
  case numel(t_psi)
    zshifteff=interp1(t_psi,zshifteff,t_th);
  case numel(t_th)
    % ok
  otherwise
    if (gdat_params.nverbose>=1)
      disp(' bad time dimension for zshift')
      disp(['it should be 1 or ' num2str(numel(t_th)) ' or ' num2str(numel(t_psi))])
      return
    end
end

% Get LIUQE times corresponding to TS times
t_eq = t_psi(iround(t_psi,t_th));

% Get PSI map
psi = psitbximas(gdat_data.shot,t_eq,'*0',psitbx_str);
% PSITBXIMAS will have removed duplicate times
i_psi = iround(psi.t,t_eq); % Mapping from psi times to TS times

% Define grid points where to evaluate PSI
if ~any(diff(zshifteff))
  % Constant vector
  pgrid = psitbxgrid('cylindrical','points',{th_r,th_z - zshifteff(1),0});
else
  % We need as many time points in psi as TS times
  psi = subs_time(psi,i_psi); % Now has ntt times
  i_psi = 1:ntt; % Trivial mapping from psi times to TS times
  % Time-varying vector
  th_z_eff = repmat(th_z(:),1,ntt) - repmat(reshape(zshifteff,1,ntt),nch,1);
  th_r_eff = repmat(th_r(:),1,ntt);
  pgrid = psitbxgrid('cylindrical','time-points',{th_r_eff,th_z_eff,0});
end

% Compute psi at TS positions
psi_ts = psitbxf2f(psi,pgrid);
psiscatvol.data = squeeze(psi_ts.x(:,i_psi));
psiscatvol.dim{1} = gdat_data.x;
psiscatvol.dim{2} = t_th;
% NOTE: we should probably not include time points where equilibrium time is far from TS time.

% Compute psi_axis at TS times
% Do not use psi_axis node because:
% - For legacy LIUQE, psi_axis is not computed by psitbx, so we could get a
% different answer and cause complex numbers computing rho_psi
% - For MATLAB LIUQE, psi_axis is only the result of asxymex whereas the
% psitbx adds some Newton iterations so again complex numbers are possible
psi_norm = psitbxp2p(psi,'01');
psi_max.data = psi_norm.psimag(i_psi);
psi_max.dim = {t_th};