gdat_tcv.m

function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(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
%
% Inputs:
%
%    no inputs: return default parameters in a structure form in gdat_params
%    shot: shot number
%    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','TCV' (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...))
% 
%    [a1,a2]=gdat;
%    a2.data_request = 'Ip';
%    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,'ip'); % standard call
%    a6=gdat(48836,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (note all lowercase in output)

%
% Comments for local developer:
% This 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 = 'tcv';

gdat_params.machine=default_machine;
gdat_params.doplot = 0;
gdat_params.liuqe = 11;
gdat_params.nverbose = 1;

% construct list of keywords from global set of keywords and specific TCV 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 TCV specific to all:
if ~isempty(data_request_names.tcv)
  tcv_names = fieldnames(data_request_names.tcv);
  for i=1:length(tcv_names)
    data_request_names.all.(tcv_names{i}) = data_request_names.tcv.(tcv_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 && ischar(data_request); data_request = lower(data_request); end

gdat_data.gdat_request = data_request_names_all; % so if return early gives list of possible request names
gdat_data.gdat_params.help = tcv_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)
    % means mdsopen(shot) already performed
    shot = mdsipmex(2,'$SHOT');
    gdat_data.shot = shot;
    do_mdsopen_mdsclose = 0;
    if ischar(shot) || isempty(shot)
      if gdat_params.nverbose>=1
        if isstruct(data_request) && isfield(data_request,'data_request')
          warning(['shot cannot be opened with ' data_request.data_request]);
        elseif ischar(data_request)
          warning(['shot cannot be opened with ' data_request]);
        else
          warning(['shot cannot be opened']);
        end
      end
      return
    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.data_request = lower(data_request.data_request);
    data_request_eff = 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})
        % enforce lower case for any character driven input
        if ischar(varargin_eff{i+1}) && ~strcmp('path',varargin_eff{i})
          gdat_params.(lower(varargin_eff{i})) = lower(varargin_eff{i+1});
        else
          gdat_params.(lower(varargin_eff{i})) = varargin_eff{i+1};
        end
      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
data_request_eff = gdat_params.data_request; % in case was defined in pairs

% if it is a request_keyword copy it:
if ischar(data_request_eff) || length(data_request_eff)==1
  ij=strmatch(data_request_eff,data_request_names_all,'exact');
else
  ij=[];
end
if ~isempty(ij); 
  gdat_data.gdat_request = data_request_names_all{ij};
  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;

% 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;
error_status = 6; % at least reached this level

liuqe_version = 1;
if isfield(gdat_data.gdat_params,'liuqe') && ~isempty(gdat_data.gdat_params.liuqe)
  liuqe_version = gdat_data.gdat_params.liuqe;
else
  gdat_data.gdat_params.liuqe = liuqe_version;
end
liuqe_matlab = 1; % now default should be matlab liuqe nodes
if liuqe_version<0 || (liuqe_version > 10 && liuqe_version < 20)
  liuqe_matlab = 0;
end
liuqe_version_eff = mod(liuqe_version,10);
substr_liuqe = '';
substr_liuqe_tcv_eq = '';
if liuqe_version_eff==2 || liuqe_version_eff==3 || (liuqe_matlab==1 && liuqe_version_eff==1)
  substr_liuqe = ['_' num2str(liuqe_version_eff)];
end
if liuqe_version_eff==2 || liuqe_version_eff==3
  substr_liuqe_tcv_eq = num2str(liuqe_version_eff);
end
  
% special treatment for model shot=-1 or preparation shot >=100'000
begstr = '';
if shot==-1 || (shot>=100000 && shot < 200000) || liuqe_version==-1
  % requires FBTE
  liuqe_version_eff = -1;
  begstr = 'tcv_eq( "';
  substr_liuqe = '", "FBTE" )';
end

% should replace all above by just psitbx_str...
liuqe_matlab = 1;
switch liuqe_version
 case {-1}, liuqe_ext=''; psitbx_str='FBTE'; liuqe_matlab = 0;
 case {1,21}, liuqe_ext=''; psitbx_str='LIUQE.M';
 case {11}, liuqe_ext=''; psitbx_str='LIUQE';liuqe_matlab = 0;
 case {2, 3, 22, 23}, liuqe_ext=['_' num2str(mod(liuqe_version,10))]; psitbx_str=['LIUQE.M' num2str(mod(liuqe_version,10))];
 case {12,13}, liuqe_ext=['_' num2str(mod(liuqe_version,10))]; psitbx_str=['LIUQE' num2str(mod(liuqe_version,10))];liuqe_matlab = 0;
 otherwise, error(['Unknown LIUQE version, liuqe = ' num2str(liuqe_version)]);
end

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

% Specifications on how to get the data provided in tcv_requests_mapping
mapping_for_tcv = tcv_requests_mapping(data_request_eff,shot);
gdat_data.label = mapping_for_tcv.label;
ishot=NaN;
if do_mdsopen_mdsclose
  % mdsdefaultserver tcv1.epfl.ch; % should be in tcv general path, but set-it in the meantime...
%%%  if liuqe_version_eff==-1
  if shot==-1
    ishot = mdsopen('pcs', shot);
  else
    if length(data_request_eff)>7 && strcmp(data_request_eff(1:6),'\rtc::')
      ishot = mdsopen('rtc',shot); % sub-tree
      data_request_eff = data_request_eff(7:end);
      mapping_for_tcv.expression = data_request_eff;
    else
      ishot = mdsopen(shot); % if ishot equal to shot, then mdsclose at the end
    end
  end
  if isempty(ishot) || ishot~=shot
    if gdat_params.nverbose>=1; warning(['cannot open shot= ' num2str(shot)]); end
    return
  end
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1st treat the simplest method: "tdi" (and tdiliuqe)
if strcmp(mapping_for_tcv.method(1:3),'tdi')
  % need to treat liuqe2, model, etc from options....
  substr_tdi = '';
  if liuqe_matlab==0 && strcmp(mapping_for_tcv.method,'tdiliuqe'); substr_tdi = substr_liuqe; end
  if iscell(mapping_for_tcv.expression)
    if length(mapping_for_tcv.expression)>0
      % series of arguments for tdi given in each cell
      if liuqe_matlab==1
        ij = findstr(mapping_for_tcv.expression{1},'equil_');
        if ~isempty(ij)
          mapping_for_tcv.expression{1}(ij+5:ij+6) = substr_liuqe;
        end
      end
      eval_expr = ['tdi(''' mapping_for_tcv.expression{1} ''''];
      for i=2:length(mapping_for_tcv.expression)
        eval_expr = [eval_expr ',''' mapping_for_tcv.expression{i} ''''];
      end
      eval_expr = [eval_expr ');'];
      aatmp = eval(eval_expr);
    else
      % empty or wrong expression
      error_status=701;
      return
    end
  else
    do_liuqem2liuqef = 0;
    if liuqe_version_eff==-1
      mapping_for_tcv_expression_eff = mapping_for_tcv.expression;
      if length(mapping_for_tcv.expression)>8 && strcmp(lower(mapping_for_tcv.expression(1:8)),'\results')
        mapping_for_tcv_expression_eff = mapping_for_tcv.expression(11:end);
      elseif findstr('tcv_eq',lower(mapping_for_tcv.expression))
        ij = regexpi(mapping_for_tcv.expression,'''''');
        if length(ij)<2
          disp(['expected more double quotes in mapping_for_tcv.expression = ' mapping_for_tcv.expression]);
          return
        else
          mapping_for_tcv_expression_eff = mapping_for_tcv.expression(ij(1)+2:ij(2)-1);
        end
      end
      eval_expr = ['tdi(''' begstr mapping_for_tcv_expression_eff substr_liuqe ''');']
    else
      if liuqe_matlab==1
        ij = findstr(mapping_for_tcv.expression,'equil_');
        if ~isempty(ij)
          mapping_for_tcv.expression(ij+5:ij+6) = substr_liuqe;
        end
        ij = regexpi(mapping_for_tcv.expression,'LIUQE.*\.M','once');
        if ~isempty(ij)
          mapping_for_tcv.expression = [mapping_for_tcv.expression(1:ij+6) substr_liuqe_tcv_eq mapping_for_tcv.expression(ij+7:end)];
        end
      else
        ij = regexpi(mapping_for_tcv.expression,'LIUQE.*\.M','once');
        if ~isempty(ij)
          mapping_for_tcv.expression = [mapping_for_tcv.expression(1:ij+4) substr_liuqe_tcv_eq mapping_for_tcv.expression(ij+7:end)];
          do_liuqem2liuqef = 1;
        end
      end
      eval_expr = ['tdi(''' mapping_for_tcv.expression ''');']
    end
    % special cases for matching liuqe fortran and matlab
    if liuqe_matlab == 0 && do_liuqem2liuqef==1
      % assume tcv_eq(''keyword'',''LIUQE..'')
      liuqe_fortran_matlab = liuqefortran2liuqematlab;
      ij = regexpi(eval_expr,'''''');
      if numel(ij)>=2
        liuqe_keyword = eval_expr(ij(1)+2:ij(2)-1);
      end
      ij_row=strmatch(liuqe_keyword,liuqe_fortran_matlab(:,2),'exact');
      if ~isempty(ij_row)
        eval_expr = [eval_expr(1:ij(1)+1) liuqe_fortran_matlab{ij_row,1} eval_expr(ij(1)+2+length(liuqe_keyword):end)];
      end
    end
    aatmp=eval(eval_expr);
  end
  if isempty(aatmp.data) || (isempty(aatmp.dim) && ischar(aatmp.data) && ~isempty(strfind(lower(aatmp.data),'no data')))% || ischar(aatmp.data) (to add?)
    if (gdat_params.nverbose>=1); warning(['problems loading data for ' eval_expr ' for data_request= ' data_request_eff]); end
    if (gdat_params.nverbose>=3); disp('check .gdat_request list'); end
    return
  end
  gdat_data.data = aatmp.data;
  gdat_data.dim = aatmp.dim;
  nbdims = length(gdat_data.dim);
  if mapping_for_tcv.timedim==-1; 
    % try to find time dim from units
    idim_non1 = []; len_dim = [];
    for i=1:length(aatmp.dimunits)
      if strcmp(aatmp.dimunits{i},'s'); mapping_for_tcv.timedim = i; end
      if length(aatmp.dim{i})>1
        idim_non1(end+1) = i;
        len_dim(end+1) = length(aatmp.dim{i});
      end
    end
    if length(idim_non1)==1
      % only one dim non 1, assume it is time
      mapping_for_tcv.timedim = idim_non1(1);
    else
      [aamax,iaamax]=max(len_dim);
      if aamax./min(len_dim)>100
        mapping_for_tcv.timedim = idim_non1(iaamax);
      end
    end
    if mapping_for_tcv.timedim==-1
      % assume last one except if of length 1
      mapping_for_tcv.timedim = nbdims;
      if (nbdims>1 && size(gdat_data.data,nbdims)==1); mapping_for_tcv.timedim = nbdims-1; end
    end
  end
  dim_nontim = setdiff([1:nbdims],mapping_for_tcv.timedim);
  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
  if mapping_for_tcv.timedim>0; gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim}; end
  gdat_data.units = aatmp.units;
  gdat_data.dimunits = aatmp.dimunits;
  if mapping_for_tcv.gdat_timedim>0 && mapping_for_tcv.gdat_timedim ~= mapping_for_tcv.timedim
    % shift timedim to gdat_timedim data(i,j,...itime,k,...) -> data(i,inewtime,j,...,k,...)
    % note that this means that gdat_data.x and gdat_data.t are same and correct, 
    % only .data, .dim and .dimunits need to be changed
    iprev=[1:nbdims];
    ij=find(dim_nontim>mapping_for_tcv.gdat_timedim-1);
    inew=[1:mapping_for_tcv.gdat_timedim-1 mapping_for_tcv.timedim dim_nontim(ij)];
    data_sizes = size(aatmp.data);
    gdat_data.data = NaN*ones(data_sizes(inew));
    abcol=ones(1,nbdims)*double(':'); abcomma=ones(1,nbdims)*double(',');
    dimstr_prev=['(' repmat(':,',1,mapping_for_tcv.timedim-1) 'it,' ...
                 repmat(':,',1,nbdims-mapping_for_tcv.timedim-1) ':)'];
    dimstr_new=['(' repmat(':,',1,mapping_for_tcv.gdat_timedim-1) 'it,' ...
                repmat(':,',1,nbdims-mapping_for_tcv.gdat_timedim-1) ':)'];
    % eval gdat_data.data(;,:,...,it,...) = aatmp.data(:,:,:,it,...);
    for it=1:size(aatmp.data,mapping_for_tcv.timedim)
      shift_eval = ['gdat_data.data' dimstr_new ' = aatmp.data'  dimstr_prev ';'];
      eval(shift_eval);
    end
    gdat_data.dim = aatmp.dim(inew);
    gdat_data.dimunits = aatmp.dimunits(inew);
  else
    mapping_for_tcv.gdat_timedim = mapping_for_tcv.timedim;
  end
  gdat_data.data_fullpath=[mapping_for_tcv.expression];
  gdat_data.help = aatmp.help;

  % end of method "tdi"

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_tcv.method,'expression')
  % 2nd: 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_tcv.expression ';'];
  eval([mapping_for_tcv.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_tcv.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_tcv.timedim==-1; 
      mapping_for_tcv.timedim = nbdims;
      if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_tcv.timedim = nbdims-1; end
    end
    dim_nontim = setdiff([1:nbdims],mapping_for_tcv.timedim);
    ijt=find(strcmp(tmp_fieldnames,'t')==1);
    if isempty(ijt)
      gdat_data.t = gdat_data.dim{mapping_for_tcv.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_tcv.expression;
    gdat_data.help = gdat_tmp.help;
  end
  % end of method "expression"

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_tcv.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 {'a_minor','rgeom','r_geom','a_minor_rho','r_geom_rho','rgeom_rho'}
    if strcmp(data_request_eff,'r_geom'); data_request_eff = 'rgeom'; end 
    if strcmp(data_request_eff,'r_geom_rho'); data_request_eff = 'rgeom_rho'; end 
    % compute average minor or major radius (on z=zaxis normally)
    if liuqe_matlab==0
      nodenameeff=['tcv_eq(''r_max_psi'',''LIUQE' substr_liuqe_tcv_eq ''')'];
      rmaxpsi=tdi(nodenameeff);
      ijnan = find(isnan(rmaxpsi.data));
      if isempty(rmaxpsi.data) || isempty(rmaxpsi.dim) || ischar(rmaxpsi.data) || ...
            ( ~isempty(ijnan) && prod(size(ijnan))==prod(size(rmaxpsi.data)) )
        if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
        if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
        return
      end   
      nodenameeff2=['tcv_eq(''r_min_psi'',''LIUQE' substr_liuqe_tcv_eq ''')'];
      rminpsi=tdi(nodenameeff2);
      ijnan = find(isnan(rminpsi.data));
      if isempty(rminpsi.data) || isempty(rminpsi.dim) || ...
            ( ~isempty(ijnan) && prod(size(ijnan))==prod(size(rminpsi.data)) )
        if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff2 ' for data_request= ' data_request_eff]); end
        if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
        return
      end
      ij=find(rmaxpsi.data<0.5 | rmaxpsi.data>1.2); % points outside TCV vessel
      if ~isempty(ij); rmaxpsi.data(ij)=NaN; end
      ij=find(rminpsi.data<0.5 | rminpsi.data>1.2);
      if ~isempty(ij); rminpsi.data(ij)=NaN; end
      if strcmp(data_request_eff,'a_minor')
        gdat_data.data=0.5.*(rmaxpsi.data(end,:) - rminpsi.data(end,:));
        gdat_data.data_fullpath=[nodenameeff ' - ' nodenameeff2 ' /2'];
      elseif strcmp(data_request_eff,'a_minor_rho')
        gdat_data.data=0.5.*(rmaxpsi.data(:,:) - rminpsi.data(:,:));
        gdat_data.data_fullpath=[nodenameeff ' - ' nodenameeff2 ' /2'];
      elseif strcmp(data_request_eff,'rgeom')
        gdat_data.data=0.5.*(rmaxpsi.data(end,:) + rminpsi.data(end,:));
        gdat_data.data_fullpath=[nodenameeff ' + ' nodenameeff2 ' /2'];
      elseif strcmp(data_request_eff,'rgeom_rho')
        gdat_data.data=0.5.*(rmaxpsi.data(:,:) + rminpsi.data(:,:));
        gdat_data.data_fullpath=[nodenameeff ' + ' nodenameeff2 ' /2'];
      else
        if (gdat_params.nverbose>=1); warning(['should not be in this case with data_request_eff = ' data_request_eff]); end
        return
      end
      if strcmp(data_request_eff(end-3:end),'_rho')
        gdat_data.dim = rmaxpsi.dim;
        gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
        gdat_data.x = gdat_data.dim{1};
        gdat_data.dimunits = rmaxpsi.dimunits(2);
      else
        gdat_data.dim = rmaxpsi.dim(2);
        gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
        gdat_data.dimunits = rmaxpsi.dimunits(2);
      end
      if any(strcmp(fieldnames(rmaxpsi),'units'))
        gdat_data.units = rmaxpsi.units;
      end
    else
      if isempty(substr_liuqe); substr_liuqe = '_1'; end
      if strcmp(data_request_eff,'a_minor') % to update when ready, still buggy, this should be rho,t
        nodenameeff=['tcv_eq(''a_minor_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
      elseif strcmp(data_request_eff,'a_minor_rho')
        nodenameeff=['tcv_eq(''r_out_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
        nodenameeff2=['tcv_eq(''r_in_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
        nodenameeff=['tcv_eq(''r_out'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
        nodenameeff2=['tcv_eq(''r_in'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
      elseif strcmp(data_request_eff,'rgeom')
        nodenameeff=['tcv_eq(''r_geom_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
      elseif strcmp(data_request_eff,'rgeom_rho')
        nodenameeff=['tcv_eq(''r_out_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
        nodenameeff2=['tcv_eq(''r_in_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
        nodenameeff=['tcv_eq(''r_out'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
        nodenameeff2=['tcv_eq(''r_in'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
      else
        if (gdat_params.nverbose>=1); warning(['should not be in this case with data_request_eff = ' data_request_eff]); end
        return
      end
      aatmp = tdi(nodenameeff);
      if strcmp(data_request_eff(end-3:end),'_rho')
        aatmp2 = tdi(nodenameeff2);
        if strcmp(data_request_eff,'a_minor_rho')
          gdat_data.data = 0.5.*(aatmp.data-aatmp2.data);
        elseif strcmp(data_request_eff,'rgeom_rho')
          gdat_data.data = 0.5.*(aatmp.data+aatmp2.data);
        end
        gdat_data.dim = aatmp.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes
        gdat_data.x = gdat_data.dim{1};
        gdat_data.data_fullpath=[nodenameeff ' and ' nodenameeff2];
      else
        gdat_data.data = aatmp.data(end,:)';
        gdat_data.dim = aatmp.dim(2); 
        aatmp.dimunits = aatmp.dimunits(2); 
        gdat_data.data_fullpath=[nodenameeff];
      end
      gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
      gdat_data.units = aatmp.units(end);
      gdat_data.dimunits = aatmp.dimunits;
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'zgeom','z_geom'}
    % compute average minor or major radius (on z=zaxis normally)
    %
    if liuqe_matlab==0
      nodenameeff=['tcv_eq(''z_contour'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
    else
      if isempty(substr_liuqe); substr_liuqe = '_1'; end
      nodenameeff=['tcv_eq(''z_edge'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
    end
    zcontour=tdi(nodenameeff);
    if isempty(zcontour.data) || isempty(zcontour.dim)  % || ischar(zcontour.data) (to add?)
      if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
      if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
      return
    end   
    if strcmp(data_request_eff,'zgeom')
      gdat_data.data=0.5.*(max(zcontour.data,[],1) + min(zcontour.data,[],1));
      gdat_data.data_fullpath=['(max+min)/2 of ' nodenameeff];
      gdat_data.dim{1} = zcontour.dim{2};
      gdat_data.dimunits{1} = zcontour.dimunits{2};
    else
      if (gdat_params.nverbose>=1); warning(['should not be in this case with data_request_eff = ' data_request_eff]); end
      return
    end
    gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
    if any(strcmp(fieldnames(zcontour),'units'))
      gdat_data.units = zcontour.units;
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'b0'}
    % B0 at R0=0.88
    R0EXP=0.88;
    if liuqe_version_eff==-1
      nodenameeff = 'tcv_eq("BZERO","FBTE")';
      tracetdi=tdi(nodenameeff);
      gdat_data.data = tracetdi.data;
    else
      nodenameeff=['\magnetics::iphi'];
      tracetdi=tdi(nodenameeff);
      gdat_data.data=192.E-07 * 0.996 *tracetdi.data/R0EXP;
    end
    if isempty(tracetdi.data) || isempty(tracetdi.dim) % || ischar(tracetdi.data) (to add?)
      if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
      return
    end
    gdat_data.data_fullpath=[nodenameeff];
    gdat_data.dim = tracetdi.dim;    
    gdat_data.t = gdat_data.dim{1};
    if any(strcmp(fieldnames(tracetdi),'units'))
      gdat_data.units = tracetdi.units;
    end
    gdat_data.dimunits = tracetdi.dimunits;
    gdat_data.request_description = ['vacuum magnetic field at R0=' num2str(R0EXP) 'm; COCOS=17'];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'betan'}
    % 100*beta / |Ip[MA] * B0[T]| * a[m]
    % get B0 from gdat_tcv, without re-opening the shot and using the same parameters except data_request
    % easily done thanks to structure call for options
    params_eff = gdat_data.gdat_params;
    params_eff.data_request='b0'; 
    b0=gdat_tcv([],params_eff); % note: no need to set .doplot=0 since gdat_tcv does not call gdat_plot in any case
    params_eff.data_request='ip';
    ip=gdat_tcv([],params_eff);
    params_eff.data_request='beta';
    beta=gdat_tcv([],params_eff);
    params_eff.data_request='a_minor';
    a_minor=gdat_tcv([],params_eff);
    % use beta as time base
    if isempty(b0.data) || isempty(b0.dim) || isempty(ip.data) || isempty(ip.dim) || isempty(a_minor.data) || isempty(a_minor.dim) || isempty(beta.data) || isempty(beta.dim)
      if (gdat_params.nverbose>=1); warning(['problems loading data for data_request= ' data_request_eff]); end
      return
    end
    gdat_data.dim = beta.dim;
    gdat_data.t = beta.dim{1};
    gdat_data.data = beta.data;
    ij=find(isfinite(ip.data));
    ip_t = interp1(ip.dim{1}(ij),ip.data(ij),gdat_data.t);
    ij=find(isfinite(b0.data));
    b0_t = interp1(b0.dim{1}(ij),b0.data(ij),gdat_data.t);
    ij=find(isfinite(a_minor.data));
    a_minor_t = interp1(a_minor.dim{1}(ij),a_minor.data(ij),gdat_data.t);
    gdat_data.data = 100.*beta.data ./ abs(ip_t).*1.e6 .* abs(b0_t) .* a_minor_t;
    gdat_data.data_fullpath='100*beta/ip*1e6*b0*a_minor, each from gdat_tcv';
    gdat_data.units = '';
    gdat_data.dimunits{1} = 's';
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'cxrs','cxrs_rho'}
    %not yet finished, just started
    % load typical data from cxrs, Ti, ni, vtori and vpoli (if available), as well as zeff from cxrs
    % if 'fit' option is added: 'fit',1, then the fitted profiles are returned
    % 
    % sub_nodes names from CXRS_get_profiles function, lower case when used in gdat
    sub_nodes = {'Ti','vTor','vPol','nC','Zeff'}; % first node is also copied into data, choose "default' one
    sub_nodes_out = lower({'Ti','vTor','vPol','ni','Zeff'}); % 
    sub_nodes_units = {'eV','m/s','m/s','m^{-3}',''}; % first node is also copied into data, choose "default' one
    % use A. Karpushov routine to get profiles and then copy the data or the fitted profiles
    aa=CXRS_get_profiles; cxrs_params = aa.param;
    cxrs_params.k_plot=0; cxrs_params.k_debug=0;
    % add params from gdat call
    params_eff = gdat_data.gdat_params;
    if isfield(params_eff,'cxrs_plot') && params_eff.cxrs_plot>0
      cxrs_plot = params_eff.cxrs_plot;
    else
      cxrs_plot = 0;
    end
    gdat_data.gdat_params.cxrs_plot = cxrs_plot;
    if isfield(params_eff,'source') && ~isempty(params_eff.source)
      source = params_eff.source;
    else
      source = [1 2 3];
    end
    gdat_data.gdat_params.source = source;
    if isfield(params_eff,'time_interval') && ~isempty(params_eff.time_interval) && length(params_eff.time_interval)>=2
      time_interval = params_eff.time_interval;
      cxrs_plot=1;
    else
      time_interval = [];
    end
    gdat_data.gdat_params.time_interval = time_interval;
    gdat_data.gdat_params.cxrs_plot = cxrs_plot;
    fit_tension_default = -100.;
    if isfield(params_eff,'fit_tension')
      fit_tension = params_eff.fit_tension;
    else
      fit_tension = fit_tension_default;
    end
    if ~isstruct(fit_tension)
      fit_tension_eff.ti = fit_tension;
      fit_tension_eff.vi = fit_tension;
      fit_tension_eff.ni = fit_tension;
      fit_tension_eff.zeff = fit_tension;
      fit_tension = fit_tension_eff;
    else
      if ~isfield(fit_tension,'ti'); fit_tension.ti = fit_tension_default; end
      if ~isfield(fit_tension,'vi'); fit_tension.vi = fit_tension_default; end
      if ~isfield(fit_tension,'ni') && ~isfield(fit_tension,'nc'); fit_tension.ni = fit_tension_default; end
      if ~isfield(fit_tension,'zeff'); fit_tension.zeff = fit_tension_default; end
    end
    gdat_data.gdat_params.fit_tension = fit_tension;
    cxrs_params.prof.Ti.taus = fit_tension.ti;
    cxrs_params.prof.vi.taus = fit_tension.vi;
    cxrs_params.prof.nc.taus = fit_tension.ni;
    cxrs_params.prof.zeff.taus = fit_tension.zeff; 
    cxrs_params.k_plot = cxrs_plot;
    cxrs_profiles = CXRS_get_profiles(shot,source,time_interval,cxrs_params);
    inb_times = length(cxrs_profiles.Times);
    gdat_data.cxrs_params = cxrs_profiles.param;
    if isempty(cxrs_profiles.Times) || ~isfield(cxrs_profiles,'proffit')
      if (gdat_params.nverbose>=1); warning(['problems loading data with CXRS_get_profiles for data_request= ' data_request_eff]); end
      for i=1:length(sub_nodes)
        sub_eff_out = sub_nodes_out{i};
        gdat_data.(sub_eff_out).fit.data = [];
        gdat_data.(sub_eff_out).fit.rho = [];
        gdat_data.(sub_eff_out).fit.error_bar = [];
        gdat_data.(sub_eff_out).raw.data = [];
        gdat_data.(sub_eff_out).raw.rho = [];
        gdat_data.(sub_eff_out).raw.error_bar = [];
        gdat_data.(sub_eff_out).raw.error_bar_rho = [];
        gdat_data.(sub_eff_out).raw.cxrs_system = [];
        gdat_data.time_interval = [];
        gdat_data.(sub_eff_out).units = sub_nodes_units{i};
        if i==1
          gdat_data.data = gdat_data.(sub_eff_out).fit.data;
          gdat_data.x = gdat_data.(sub_eff_out).fit.rho;
          gdat_data.error_bar = gdat_data.(sub_eff_out).fit.error_bar;
          gdat_data.units = gdat_data.(sub_eff_out).units;
        end
      end
      return
    end
    inb_channels =120; % need to change if gets bigger!!! but easier to prefill with NaNs and use the "use" part
    for i=1:length(sub_nodes)
      sub_eff = sub_nodes{i};
      sub_eff_out = sub_nodes_out{i};
      % fits
      if isfield(cxrs_profiles.proffit,sub_eff)
        gdat_data.(sub_eff_out).fit.data = cxrs_profiles.proffit.(sub_eff);
        gdat_data.(sub_eff_out).fit.rho = cxrs_profiles.proffit.([sub_eff '_rho']);
        gdat_data.(sub_eff_out).fit.error_bar = cxrs_profiles.proffit.(['d' sub_eff]);
      else
        gdat_data.(sub_eff_out).fit.data = [];
        gdat_data.(sub_eff_out).fit.rho = [];
        gdat_data.(sub_eff_out).fit.error_bar = [];
      end
      % raw data (use all data so keep same size)
      gdat_data.(sub_eff_out).raw.data = NaN * ones(inb_channels,inb_times);
      gdat_data.(sub_eff_out).raw.rho = NaN * ones(inb_channels,inb_times);
      gdat_data.(sub_eff_out).raw.error_bar = NaN * ones(inb_channels,inb_times);
      gdat_data.(sub_eff_out).raw.error_bar_rho = NaN * ones(inb_channels,inb_times);
      gdat_data.(sub_eff_out).raw.cxrs_system = NaN * ones(inb_channels,inb_times);
      gdat_data.time_interval = [];
      for it=1:inb_times
        if isfield(cxrs_profiles,sub_eff)
          nb_raw_points = length(cxrs_profiles.(sub_eff){it}.use.y);
          gdat_data.(sub_eff_out).raw.data(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.y;
          gdat_data.(sub_eff_out).raw.rho(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.x;
          gdat_data.(sub_eff_out).raw.error_bar(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.dy;
          gdat_data.(sub_eff_out).raw.error_bar_rho(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.dx;
          gdat_data.(sub_eff_out).raw.cxrs_system(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.sys;
          gdat_data.time_interval{it} = cxrs_profiles.(sub_eff){it}.t_lim;
        end          
      end
      gdat_data.(sub_eff_out).units = sub_nodes_units{i};
      if i==1
        gdat_data.data = gdat_data.(sub_eff_out).fit.data;
        gdat_data.x = gdat_data.(sub_eff_out).fit.rho;
        gdat_data.error_bar = gdat_data.(sub_eff_out).fit.error_bar;
        gdat_data.units = gdat_data.(sub_eff_out).units;
      end
    end
    gdat_data.t = cxrs_profiles.proffit.time;
    gdat_data.dim = {gdat_data.x; gdat_data.t};
    if isempty(time_interval)
      gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[1 2 3],[],cxrs_params); % with cxrs_params'];
    else
      gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[1 2 3],[' num2str(time_interval) '],cxrs_params); % with cxrs_params'];
    end
    gdat_data.dimunits{1} = '';
    gdat_data.dimunits{2} = 's';
    % add grids_1d to have rhotor, etc if cxrs_rho was asked for
    if strcmp(data_request_eff,'cxrs_rho')
      gdat_data = get_grids_1d(gdat_data,2,1,gdat_params.nverbose);
    else
      gdat_data = get_grids_1d(gdat_data,2,0,gdat_params.nverbose);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'eqdsk'}
    %
    time=1.; % default time
    if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time)
      time = gdat_data.gdat_params.time;
    else
      gdat_data.gdat_params.time = time;
      if (gdat_params.nverbose>=3); disp(['"time" is expected as an option, choose default time = ' num2str(time)]); end
    end
    gdat_data.gdat_params.time = time;
    gdat_data.t = time;
    zshift = 0.;
    if isfield(gdat_data.gdat_params,'zshift') && ~isempty(gdat_data.gdat_params.zshift)
      zshift = gdat_data.gdat_params.zshift;
    else
      gdat_data.gdat_params.zshift = zshift;
    end
    gdat_data.gdat_params.zshift = zshift;
    nrz_out = [129,129];
    if isfield(gdat_data.gdat_params,'nrz_out') && ~isempty(gdat_data.gdat_params.nrz_out)
      nrz_out = gdat_data.gdat_params.nrz_out;
    else
      gdat_data.gdat_params.nrz_out = nrz_out;
    end
    gdat_data.gdat_params.nrz_out = nrz_out;
    for itime=1:length(time)
      time_eff = time(itime);
      % use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2
      [fnames_readresults]=read_results_for_chease(shot,time_eff,liuqe_version,3,[],[],[],zshift,0,1,nrz_out);
      cocos_read_results_for_chease = 17;
      if isempty(fnames_readresults)
        if gdat_params.nverbose>=1;
          warning(['could not create eqdsk file from read_results_for_chease with data_request= ' data_request_eff]);
        end
        return
      end
      eqdskval=read_eqdsk(fnames_readresults{4},cocos_read_results_for_chease,0,[],[],1); % read_results now has relevant cocos LIUQE=17
      for i=1:length(fnames_readresults)
        unix(['rm ' fnames_readresults{i}]);
      end
      % transform to cocos=2 since read_results originally assumed it was cocos=2
      cocos_in = 2;
      [eqdsk_cocos_in, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdskval,[cocos_read_results_for_chease cocos_in]);
      fnamefull = fullfile(['/tmp/' getenv('USER')],['EQDSK_' num2str(shot) 't' num2str(time_eff,'%.4f')]);
      % We still write COCOS=2 case, since closer to standard (in /tmp)
      write_eqdsk(fnamefull,eqdsk_cocos_in,cocos_in,[],[],[],1);
      % Now gdat_tcv should return the convention from LIUQE which is COCOS=17, except if specified in option
      % create standard filename name from shot, time_eff (cocos will be added by write_eqdsk)
      cocos_out = 17;
      if isfield(gdat_data.gdat_params,'cocos') && ~isempty(gdat_data.gdat_params.cocos)
        cocos_out = gdat_data.gdat_params.cocos;
      else
        gdat_data.gdat_params.cocos = cocos_out;
      end
      [eqdsk_cocosout, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdsk_cocos_in,[cocos_in cocos_out]);
      % for several times, use array of structure for eqdsks, 
      % cannot use it for psi(R,Z) in .data and .x since R, Z might be different at different times,
      % so project psi(R,Z) on Rmesh, Zmesh of 1st time
      if length(time) > 1
        gdat_data.eqdsk{itime} = write_eqdsk(fnamefull,eqdsk_cocosout,cocos_out);
        if itime==1
          gdat_data.data(:,:,itime) = gdat_data.eqdsk{itime}.psi;
          gdat_data.dim{1} = gdat_data.eqdsk{itime}.rmesh;
          gdat_data.dim{2} = gdat_data.eqdsk{itime}.zmesh;
        else
	  xx=repmat(reshape(gdat_data.dim{1},length(gdat_data.dim{1}),1),1,size(gdat_data.eqdsk{itime}.psi,2));
	  yy=repmat(reshape(gdat_data.dim{2},1,length(gdat_data.dim{2})),size(gdat_data.eqdsk{itime}.psi,1),1);
	  aa = interpos2Dcartesian(gdat_data.eqdsk{itime}.rmesh,gdat_data.eqdsk{itime}.zmesh ...
	  ,gdat_data.eqdsk{itime}.psi,xx,yy,-1,-1);
          gdat_data.data(:,:,itime) = aa;
        end
      else
        gdat_data.eqdsk = write_eqdsk(fnamefull,eqdsk_cocosout,cocos_out);
        gdat_data.data = gdat_data.eqdsk.psi;
        gdat_data.dim{1} = gdat_data.eqdsk.rmesh;
        gdat_data.dim{2} = gdat_data.eqdsk.zmesh;
      end
    end
    gdat_data.dim{3} = gdat_data.t;
    gdat_data.x = gdat_data.dim(1:2);
    gdat_data.data_fullpath=['psi(R,Z) and eqdsk from read_eqdsk from LIUQE' num2str(liuqe_version) ';zshift=' num2str(zshift)];
    gdat_data.units = 'T m^2';
    gdat_data.dimunits = {'m','m','s'};
    gdat_data.request_description = ['data=psi, x=(R,Z), eqdsk contains eqdsk structure with which ' ...
                    'plot_eqdsk, write_eqdsk, read_eqdsk can be used'];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'halpha'}
    channels = -1;
    if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels)
      channels = gdat_data.gdat_params.channels;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'mhd'}
    % load n=1, 2 and 3 Bdot from magnetic measurements
    if shot< 50926
      n1=tdi('abs(mhdmode("LFS",1,1))');
      n2=tdi('abs(mhdmode("LFS",2,1))');
      n3=tdi('abs(mhdmode("LFS",3,1))');
      gdat_data.data_fullpath='abs(mhdmode("LFS",n,1));% for 1, 2, 3';
    else
      if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source)
        % gdat_data.gdat_params.source;
      else
        gdat_data.gdat_params.source = '23';
      end      
      if length(gdat_data.gdat_params.source)>=2 && strcmp(gdat_data.gdat_params.source(1:2),'23')
        aaLFSz23_sect3=tdi('\atlas::DT196_MHD_001:channel_067');
        aaLFSz23_sect11=tdi('\atlas::DT196_MHD_001:channel_075');
        n1 = aaLFSz23_sect3;
        n1.data = aaLFSz23_sect3.data - aaLFSz23_sect11.data;
        n2 = aaLFSz23_sect3;
        n2.data = aaLFSz23_sect3.data + aaLFSz23_sect11.data;
        n3=n1;
        gdat_data.data_fullpath='\atlas::DT196_MHD_001:channel_067 -+ \atlas::DT196_MHD_001:channel_075 for n=1,2, LFS_sect_3/11, z=23cm';
        if strcmp(gdat_data.gdat_params.source,'23full')
          % HFS from sec 3 and 11
          aaHFSz23_sect3=tdi('\atlas::DT196_MHD_002:channel_018');
          aaHFSz23_sect11=tdi('\atlas::DT196_MHD_002:channel_022');
          gdat_data.n1_HFS=aaHFSz23_sect3;
          gdat_data.n1_HFS.data = gdat_data.n1_HFS.data - aaHFSz23_sect11.data;
          gdat_data.n2_HFS=aaHFSz23_sect3;
          gdat_data.n2_HFS.data = gdat_data.n2_HFS.data + aaHFSz23_sect11.data;
          gdat_data.data_fullpath=['\atlas::DT196_MHD_001:channel_067 -+ \atlas::DT196_MHD_001:channel_075 for n=1,2, LFS_sect_3/11, z=23cm; _HFS' ...
                    ' same for sector HFS: MHD_002:channel_018-+MHD_002:channel_022'];
        end
      elseif strcmp(gdat_data.gdat_params.source(1),'0')
        aaLFSz0_sect3=tdi('\atlas::DT196_MHD_001:channel_083');
        aaLFSz0_sect11=tdi('\atlas::DT196_MHD_001:channel_091');
        n1 = aaLFSz0_sect3;
        n1.data = aaLFSz0_sect3.data - aaLFSz0_sect11.data;
        n2 = aaLFSz0_sect3;
        n2.data = aaLFSz0_sect3.data + aaLFSz0_sect11.data;
        n3=n1;
        gdat_data.data_fullpath='\atlas::DT196_MHD_001:channel_083 -+ \atlas::DT196_MHD_001:channel_091 for n=1,2, LFS_sect_3/11, z=0cm';
        if strcmp(gdat_data.gdat_params.source,'0full')
          % sect 11 180deg from sec 3
          aaHFSz0_sect3=tdi('\atlas::DT196_MHD_002:channel_034');
          aaHFSz0_sect11=tdi('\atlas::DT196_MHD_002:channel_038');