function [gdat_data,gdat_params,error_status,varargout] = gdat_aug(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','AUG' (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
%
% Examples:
% (should add working examples for various machines (provides working shot numbers for each machine...))
%
%    [a1,a2]=gdat;
%    a2.data_request = 'Ip';
%    a3=gdat(32827,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(32827,'ip'); % standard call
%    a6=gdat(32827,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (note all lowercase in output)
%    eqd = gdat(33134,'eqdsk','equil','IDE','time',2,'extra_arg_sf2sig','''-ed'',2');
%
% to get parameters, similar special signals use: 'XXX:name' with XXX in {'param:', 'param-set:', 'area-base:', 'time-base:'}
%    gyrofreq=gdat(30594,{'ECS','P_sy2_g1','AUGD','param:gyr_freq'});
%    or
%    gyrofreq=gdat(30594,{'ECS','P_sy2_g1'},'special_signal',,'param:gyr_freq');
%
%    ip=gdat(shot,'ip'); % will show in ip.gdat_fullpath the effective signal used. In this case it is similar to ask:
%    ip=gdat(shot,{'FPC','IpiFP','AUGD'}); % or ip=gdat(shot,{'FPC','IpiFP'});
%
%    so any signal can be loaded by giving the diagnostic, signal name and experiment, etc
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Remote data access for AUG
% You need to make a "tunnel" in one unix session from the remote node using for example:
%   to create a tunnel to port 8001 to mds server lxmdsplus.aug.ipp.mpg.de
%      ssh -l ipp_username -L 8001:lxmdsplus.aug.ipp.mpg.de:8000 gate1.aug.ipp.mpg.de
%
% Then in another unix session on the remote node, in matlab and with the appropriate mds path do:
%   >> mdsconnect('localhost:8001')
%   >> mdsvalue('1+2') % should return 3 if correctly connected
%
%

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

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

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

gdat_params.machine=default_machine;
gdat_params.doplot = 0;
gdat_params.exp_name = 'AUGD';
gdat_params.nverbose = 1;

% construct list of keywords from global set of keywords and specific AUG 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 AUG specific to all:
if ~isempty(data_request_names.aug)
  aug_names = fieldnames(data_request_names.aug);
  for i=1:length(aug_names)
    data_request_names.all.(aug_names{i}) = data_request_names.aug.(aug_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_eff = data_request;
  data_request = data_request; % should not lower until see if keyword
end

gdat_data.gdat_request = data_request_names_all; % so if return early gives list of possible request names
gdat_data.gdat_params.help = aug_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=-999; % means not defined
    do_mdsopen_mdsclose = 0;
  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;
    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})
        varargin_eff{i} = lower(varargin_eff{i});
        % enforce lower case for any character driven input (only for 1st part of pair...)
% $$$         if ischar(varargin_eff{i+1})
          gdat_params.(varargin_eff{i}) = 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
% defaults for all
% default equil:
if ~isfield(gdat_params,'equil');
  gdat_params.equil = 'EQI';
  if isfield(gdat_params,'source') && (any(strcmp(lower(gdat_params.source),'ide')) || any(strcmp(lower(gdat_params.source),'idg')) ...
          || any(strcmp(lower(gdat_params.source),'ida')))
    gdat_params.equil = 'IDE';
  end
end
% $$$ if isfield(gdat_params,'source') && (any(strcmp(lower(gdat_params.source),'tre')) || any(strcmp(lower(gdat_params.source),'tra')))
% $$$   gdat_params.equil = 'TRE';
% $$$ end
extra_arg_sf2sig = '[]';
if isfield(gdat_params,'extra_arg_sf2sig') && ~isempty(gdat_params.extra_arg_sf2sig)
  extra_arg_sf2sig = gdat_params.extra_arg_sf2sig;
end
gdat_params.extra_arg_sf2sig = extra_arg_sf2sig;
%
special_signal = '';
if isfield(gdat_params,'special_signal') && ~isempty(gdat_params.special_signal)
  special_signal = gdat_params.special_signal;
end
gdat_params.special_signal = special_signal;

% if it is a request_keyword can obtain description:
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

% 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

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

% Specifications on how to get the data provided in aug_requests_mapping
mapping_for_aug = aug_requests_mapping(data_request_eff);
gdat_data.label = mapping_for_aug.label;

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1st treat the simplest method: "signal"
if strcmp(mapping_for_aug.method,'signal')
  exp_location = gdat_data.gdat_params.exp_name;
  if ~iscell(mapping_for_aug.expression)
    error_status = 1010;
    disp(['expects a cell array with at least 2 cells in expression, mapping_for_aug.expression = ' mapping_for_aug.expression]);
    if isfield(mapping_for_aug,'not_found') && mapping_for_aug.not_found
      disp(['data request not found in aug_requests_mapping, see .data_request_names_all to see all options'])
      gdat_data.data_request_names_all = data_request_names_all;
    end
    return
  elseif length(mapping_for_aug.expression)>=3 && ~isempty(mapping_for_aug.expression{3})
    exp_location = mapping_for_aug.expression{3};
  elseif length(mapping_for_aug.expression)>=2
    mapping_for_aug.expression{3} = exp_location;
  else
    error_status = 101;
    disp(['expects at least 2 cells in expression, mapping_for_aug.expression = ' mapping_for_aug.expression]);
    if isfield(mapping_for_aug,'not_found') && mapping_for_aug.not_found
      disp(['data request not found in aug_requests_mapping, see .data_request_names_all to see all options'])
      gdat_data.data_request_names_all = data_request_names_all;
    end
    return
  end
  % allow changing main source with simple signals, 'source' parameter relates to mapping_for_aug.expression{1}
  if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || ~ischar(gdat_data.gdat_params.source)
    gdat_data.gdat_params.source = mapping_for_aug.expression{1};
  else
    % special case for EQI,EQH,IDE instead of IDG when FPG is present value
    if strcmp(lower(gdat_data.gdat_params.source),'ide') && strcmp(lower(mapping_for_aug.expression{1}),'fpg')
      gdat_data.gdat_params.source = 'IDG';
    elseif strcmp(lower(gdat_data.gdat_params.source),'eqi') && strcmp(lower(mapping_for_aug.expression{1}),'fpg')
      gdat_data.gdat_params.source = 'GQI';
    elseif strcmp(lower(gdat_data.gdat_params.source),'eqh') && strcmp(lower(mapping_for_aug.expression{1}),'fpg')
      gdat_data.gdat_params.source = 'GQH';
    end
    mapping_for_aug.expression{1} = gdat_data.gdat_params.source;
  end
  gdat_data.gdat_params.exp_name = exp_location;
  % time interval
  time_interval = [];
  % extra args for rdaAUG_eff
  if length(mapping_for_aug.expression)>=4 && ~isempty(mapping_for_aug.expression{4})
    gdat_data.gdat_params.special_signal = mapping_for_aug.expression{4};
  end
  [aatmp,error_status]=rdaAUG_eff(shot,mapping_for_aug.expression{1},mapping_for_aug.expression{2},exp_location, ...
          time_interval,gdat_data.gdat_params.extra_arg_sf2sig,gdat_data.gdat_params.special_signal);
  if error_status~=0
    if gdat_params.nverbose>=3; disp(['error after rdaAUG in signal with data_request_eff= ' data_request_eff]); end
    return
  end
  if isfield(aatmp,'data'); gdat_data.data = aatmp.data; end
  if isfield(aatmp,'t'); gdat_data.t = aatmp.t; end
  if isfield(aatmp,'x'); gdat_data.x = aatmp.x; end
  if isfield(aatmp,'dim'); gdat_data.dim = aatmp.dim; end
  if isfield(aatmp,'dimunits'); gdat_data.dimunits = aatmp.dimunits; end
  gdat_data.data_fullpath=mapping_for_aug.expression;
  if isempty(aatmp.data) || (isempty(gdat_data.t) && isempty(gdat_data.x))
    return
  end
  % construct .dim if not present (means coming from sf2sig)
  if ~isfield(aatmp,'dim')
    len_t = numel(aatmp.t);
    len_x = numel(aatmp.x);
    len_data= size(aatmp.data);
    itime = find(len_data==len_t);
    nbdims_eff = sum([size(gdat_data.data)~=1]);
    aatmp.dim{itime} = aatmp.t;
    aatmp.dimunits{itime} = aatmp.time_aug.unit;
    if nbdims_eff >= 2
      % assume area-base is main 2nd coordinate (did not find a case with 2 area-base (2 coordinates defined in addition to time) to be able to test)
      i_x = setdiff([1:2],itime);
      aatmp.dim{i_x} = aatmp.x;
      if ~isempty(aatmp.area) && isfield(aatmp.area,'unit')
        aatmp.dimunits(i_x) = aatmp.area.unit;
      else
        aatmp.dimunits(i_x) = {''};
      end
      if isempty(aatmp.dimunits{i_x}); aatmp.dimunits{i_x} = ''; end
      for i=3:nbdims_eff
        aatmp.dim{i} = [1:size(aatmp.data,i)];
        aatmp.dimunits{i} = '';
      end
    end
    gdat_data.dim = aatmp.dim;
    gdat_data.dimunits = aatmp.dimunits;
  end
  mapping_for_aug_timedim_orig = mapping_for_aug.timedim;
  if mapping_for_aug.timedim<=0
    % need to guess timedim
    if prod(size(aatmp.data)) == length(aatmp.data)
      mapping_for_aug.timedim = 1;
    elseif length(size(aatmp.data))==2
      % 2 true dimensions
      if length(aatmp.t) == size(aatmp.data,1)
        mapping_for_aug.timedim = 1;
      else
        mapping_for_aug.timedim = 2;
      end
    else
      % more than 2 dimensions, find the one with same length to define timedim
      mapping_for_aug.timedim = find(size(aatmp.data)==length(aatmp.t));
      if length(mapping_for_aug.timedim); mapping_for_aug.timedim = mapping_for_aug.timedim(1); end
    end
    mapping_for_aug.gdat_timedim = mapping_for_aug.timedim;
  end
  if length(size(aatmp.data))==1 | (length(size(aatmp.data))==2 & size(aatmp.data,2)==1)
    gdat_data.dim=[{aatmp.t}];
    gdat_data.dimunits={'time [s]'};
  elseif length(size(aatmp.data))==2
    gdat_data.dim{mapping_for_aug.timedim} = gdat_data.t;
    gdat_data.dimunits{mapping_for_aug.timedim} = 'time [s]';
    gdat_data.dim{setdiff([1:2],mapping_for_aug.timedim)} = gdat_data.x;
  else
    gdat_data.dim{mapping_for_aug.timedim} = gdat_data.t;
    gdat_data.dimunits{mapping_for_aug.timedim} = 'time [s]';
    nbdims = length(size(gdat_data.data));
    % should copy dims from aatmp but skip time dm already set
    if numel(aatmp.dim) ~= nbdims
      error(['in gdat_aug: numel(aatmp.dim)(=',num2str(numel(aatmp.dim)),') ~= nbdims(=',num2str(nbdims),')']);
    else
      iotherdims = setdiff([1:nbdims],mapping_for_aug.timedim);
      % set 1st non-time dim to .x
      if ~isequal(aatmp.x,aatmp.dim{iotherdims(1)}); warning('aatmp.x was not 1st non-time dim'); end
      gdat_data.dim(iotherdims) = aatmp.dim(iotherdims);
      gdat_data.dimunits(iotherdims) = aatmp.dimunits(iotherdims);
      gdat_data.x = gdat_data.dim{iotherdims(1)};
      % often time as 1st dim, but would be more "our" standard as 2nd dim, at this stage change systematically
      if mapping_for_aug_timedim_orig <= 0 && mapping_for_aug.timedim==1
        % switch 1st two dimensions
        adim = gdat_data.dim;
        gdat_data.dim(1:2) = adim([2 1]);
        adimunits = gdat_data.dimunits;
        gdat_data.dimunits(1:2) = adimunits([2 1]);
        mapping_for_aug.timedim = 2;
        mapping_for_aug.gdat_timedim = mapping_for_aug.timedim;
        if numel(size(gdat_data.dim{1}))==2; gdat_data.dim{1} = gdat_data.dim{1}'; end
        gdat_data.x = gdat_data.dim{1};
        gdat_data.t = gdat_data.dim{mapping_for_aug.timedim};
        adata = gdat_data.data;
        gdat_data.data = permute(adata,[2 1 3:numel(size(adata))]);
      end
    end
  end
  nbdims = length(gdat_data.dim);
  gdat_data.units = aatmp.units;
  if mapping_for_aug.gdat_timedim>0 && mapping_for_aug.gdat_timedim ~= mapping_for_aug.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_aug.gdat_timedim-1);
    inew=[1:mapping_for_aug.gdat_timedim-1 mapping_for_aug.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_aug.timedim-1) 'it,' ...
                 repmat(':,',1,nbdims-mapping_for_aug.timedim-1) ':)'];
    dimstr_new=['(' repmat(':,',1,mapping_for_aug.gdat_timedim-1) 'it,' ...
                repmat(':,',1,nbdims-mapping_for_aug.gdat_timedim-1) ':)'];
    % eval gdat_data.data(;,:,...,it,...) = aatmp.data(:,:,:,it,...);
    for it=1:size(aatmp.data,mapping_for_aug.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_aug.gdat_timedim = mapping_for_aug.timedim;
  end
  % end of method "signal"

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_aug.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_aug.expression ';'];
  gdata_data_orig = gdat_data;
  fields_to_copy = {'gdat_params', 'gdat_request', 'label', 'data_fullpath', 'request_description', 'mapping_for'};
  eval([mapping_for_aug.expression ';']);
  for i=1:length(fields_to_copy)
    if isfield(gdata_data_orig,fields_to_copy{i})
      gdat_tmp.(fields_to_copy{i}) = gdata_data_orig.(fields_to_copy{i});
    end
  end
  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_aug.expression]); end
    error_status=801;
    return
  end
  tmp_fieldnames = fieldnames(gdat_tmp);
  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_aug.timedim==-1;
      mapping_for_aug.timedim = nbdims;
      if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_aug.timedim = nbdims-1; end
    end
    dim_nontim = setdiff([1:nbdims],mapping_for_aug.timedim);
    ijt=find(strcmp(tmp_fieldnames,'t')==1);
    if isempty(ijt)
      gdat_data.t = gdat_data.dim{mapping_for_aug.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_aug.expression;
  end
  % end of method "expression"

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_aug.method,'switchcase')
  switch data_request_eff % not lower(...) since data_request_eff should be lower case already at this stage
   case {'cxrs', 'cxrs_rho'}
    % if 'fit' option is added: 'fit',1, then the fitted profiles are returned which means "_rho" is effectively called
    if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit) || ~isnumeric(gdat_data.gdat_params.fit)
      if strcmp(data_request_eff,'cxrs')
        gdat_data.gdat_params.fit = 0;
      else
        gdat_data.gdat_params.fit = 1; % add fit as default if cxrs_rho requested (since equil takes most time)
      end
    elseif gdat_data.gdat_params.fit==1
      data_request_eff = 'cxrs_rho';
    end
    exp_name_eff = gdat_data.gdat_params.exp_name;
    sources_available = {'CEZ', 'CMZ', 'CUZ',    'COZ'};
    r_node_available = {'R_time', 'R', 'R_time', 'R_time'};
    z_node_available = {'z_time', 'z', 'z_time', 'z_time'};
    % scroll through list up to first available when no sources provided
    sources_available_for_scroll = {'CEZ', 'CUZ', 'COZ'};
    scroll_through_list = 0;
    if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
      gdat_data.gdat_params.source = sources_available_for_scroll(1);
      scroll_through_list = 1; % means scroll through sources until one is available
    elseif ischar(gdat_data.gdat_params.source)
      gdat_data.gdat_params.source = {gdat_data.gdat_params.source};
    end
    if length(gdat_data.gdat_params.source)==1
      if strcmp(upper(gdat_data.gdat_params.source{1}),'CEZ')
        gdat_data.gdat_params.source = {'CEZ'};
      elseif strcmp(upper(gdat_data.gdat_params.source{1}),'CMZ')
        gdat_data.gdat_params.source = {'CMZ'};
      elseif strcmp(upper(gdat_data.gdat_params.source{1}),'CUZ')
        gdat_data.gdat_params.source = {'CUZ'};
      elseif strcmp(upper(gdat_data.gdat_params.source{1}),'COZ')
        gdat_data.gdat_params.source = {'COZ'};
      elseif strcmp(upper(gdat_data.gdat_params.source{1}),'ALL')
        gdat_data.gdat_params.source = sources_available;
        scroll_through_list = 2; % means scroll through all sources and load all sources
      else
        warning(['source = ' gdat_data.gdat_params.source ' not expected with data_request= ' data_request_eff]);
        return
      end
    else
      sources_in = intersect(sources_available,upper(gdat_data.gdat_params.source));
      if length(sources_in) ~= length(gdat_data.gdat_params.source)
        disp('following sources not yet available, check with O. Sauter if need be')
        setdiff(upper(gdat_data.gdat_params.source),sources_available)
      end
      gdat_data.gdat_params.source = sources_in;
    end
    extra_arg_sf2sig_eff_string = '';
    if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
      extra_arg_sf2sig_eff_string = [',' gdat_data.gdat_params.extra_arg_sf2sig];
    end
    % set starting source
    i_count = 1;
    diag_name = gdat_data.gdat_params.source{i_count};
    sources_tried{i_count} = diag_name;
    iload = 1;
    iequil = 0;
    while iload==1
      ishotfile_ok = 1;
      i_source = strmatch(diag_name,sources_available);
      r_node = r_node_available{i_source};
      z_node = z_node_available{i_source};
      % R, Z positions of measurements
      try
        % eval(['[r_time]=sf2ab(diag_name,shot,r_node,''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
        % both for CEZ and CMZ, and.. Ti:1 is ok, otherwise introduce string above
        [r_time,err]=rdaAUG_eff(shot,diag_name,'Ti',exp_name_eff,[],extra_arg_sf2sig_eff_string,['area-base:' r_node ':1']);
      catch ME_R_time
        % assume no shotfile
        disp(getReport(ME_R_time))
        ishotfile_ok = 0;
      end
      if ishotfile_ok == 1
        gdat_data.r = r_time.data;
        inotok=find(gdat_data.r<=0);
        gdat_data.r(inotok) = NaN;
        try
          if [exist('sf2sig')==3]
            [z_time,err]=rdaAUG_eff(shot,diag_name,'z_time',exp_name_eff,[],extra_arg_sf2sig_eff_string,['area-base:' z_node ':1']);
          else
            [z_time,err]=rdaAUG_eff(shot,diag_name,'z_time',exp_name_eff,[],extra_arg_sf2sig_eff_string);
          end
          gdat_data.z = z_time.data;
          inotok=find(gdat_data.z<=0);
          gdat_data.z(inotok) = NaN;
        catch ME_R_time
          disp(getReport(ME_R_time))
          ishotfile_ok = 0;
        end
      else
        gdat_data.r = [];
        gdat_data.z = [];
      end
      if ishotfile_ok == 1
        try
          % eval(['[time]=sf2tb(diag_name,shot,''time'',''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
          if [exist('sf2sig')==3]
            [time,err] = rdaAUG_eff(shot,diag_name,'time',exp_name_eff,[],extra_arg_sf2sig_eff_string,'time-base:time:0');
          else
            time.data = r_time.t;
          end
          gdat_data.t = time.data;
        catch ME_R_time
          disp(getReport(ME_R_time))
          ishotfile_ok = 0;
        end
      else
        gdat_data.t = [];
      end
      gdat_data.dim{1} = {gdat_data.r , gdat_data.z};
      gdat_data.dimunits{1} = 'R, Z [m]';
      gdat_data.dim{2} = gdat_data.t;
      gdat_data.dimunits{2} = 't [s]';
      gdat_data.x = gdat_data.dim{1};
      % vrot
      if ishotfile_ok == 1
        try
          [a,error_status]=rdaAUG_eff(shot,diag_name,'vrot',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
          if isempty(a.data) || isempty(a.t) || error_status>0
            if gdat_params.nverbose>=3;
              a
              disp(['with data_request= ' data_request_eff])
            end
          end
          [aerr,e]=rdaAUG_eff(shot,diag_name,'err_vrot',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
        catch ME_vrot
          disp(getReport(ME_vrot))
          ishotfile_ok = 0;
        end
        gdat_data.vrot.data = a.data;
        gdat_data.vrot.error_bar = aerr.data;
      else
        gdat_data.vrot.data = [];
        gdat_data.vrot.error_bar = [];
      end
      a.name = 'vrot';
      if any(strcmp(fieldnames(a),'units')); gdat_data.vrot.units=a.units; end
      if any(strcmp(fieldnames(a),'name')); gdat_data.vrot.name=a.name; end
      gdat_data.vrot.label = 'vrot_tor';
      % Ti
      %     [a,e]=rdaAUG_eff(shot,diag_name,'Ti',exp_name_eff);
      %     [aerr,e]=rdaAUG_eff(shot,diag_name,'err_Ti',exp_name_eff);
      if ishotfile_ok == 1
        try
          [a,e]=rdaAUG_eff(shot,diag_name,'Ti_c',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
          [aerr,e]=rdaAUG_eff(shot,diag_name,'err_Ti_c',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
        catch ME_ti
          disp(getReport(ME_ti))
          ishotfile_ok = 0;
        end
        gdat_data.ti.data = a.data;
        gdat_data.ti.error_bar = aerr.data;
      else
        gdat_data.ti.data = [];
        gdat_data.ti.error_bar = [];
      end
      a.name = 'Ti_c';
      if any(strcmp(fieldnames(a),'units')); gdat_data.ti.units=a.units; end
      if any(strcmp(fieldnames(gdat_data.ti),'units')); gdat_data.units=gdat_data.ti.units; end
      if any(strcmp(fieldnames(a),'name')); gdat_data.ti.name=a.name; end
      gdat_data.ti.label = 'Ti_c';
      % main node ti a this stage
      gdat_data.data = gdat_data.ti.data;
      gdat_data.label = [gdat_data.label '/Ti'];
      gdat_data.error_bar = gdat_data.ti.error_bar;
      gdat_data.data_fullpath=[exp_name_eff '/' diag_name '/' a.name ';vrot, Ti in data, {r;z} in .x'];
      %
      if strcmp(data_request_eff,'cxrs_rho')
        % defaults
        if iequil == 0
          gdat_data.equil.data = [];
          gdat_data.psi = [];
          gdat_data.rhopolnorm = [];
          gdat_data.rhotornorm = [];
          gdat_data.rhovolnorm = [];
          % defaults for fits, so user always gets std structure
          gdat_data.fit.rhotornorm = []; % same for both ti and vrot
          gdat_data.fit.rhopolnorm = [];
          gdat_data.fit.t = [];
          gdat_data.fit.ti.data = [];
          gdat_data.fit.ti.drhotornorm = [];
          gdat_data.fit.vrot.data = [];
          gdat_data.fit.vrot.drhotornorm = [];
          gdat_data.fit.raw.rhotornorm = [];
          gdat_data.fit.raw.ti.data = [];
          gdat_data.fit.raw.vrot.data = [];
          fit_tension_default = -1;
          if isfield(gdat_data.gdat_params,'fit_tension')
            fit_tension = gdat_data.gdat_params.fit_tension;
          else
            fit_tension = fit_tension_default;
          end
          if ~isstruct(fit_tension)
            fit_tension_eff.ti = fit_tension;
            fit_tension_eff.vrot = fit_tension;
            fit_tension = fit_tension_eff;
          else
            if ~isfield(fit_tension,'ti'); fit_tension.ti = fit_tension_default; end
            if ~isfield(fit_tension,'vrot '); fit_tension.vrot = fit_tension_default; end
          end
          gdat_data.gdat_params.fit_tension = fit_tension;
          if isfield(gdat_data.gdat_params,'fit_nb_rho_points')
            fit_nb_rho_points = gdat_data.gdat_params.fit_nb_rho_points;
          else
            fit_nb_rho_points = 201;
          end
          gdat_data.gdat_params.fit_nb_rho_points = fit_nb_rho_points;
        end
        if ishotfile_ok == 1 && iequil == 0
          params_equil = gdat_data.gdat_params;
          params_equil.data_request = 'equil';
          [equil,params_equil,error_status] = gdat_aug(shot,params_equil);
          if error_status>0
            if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end
            return
          end
          iequil = 1;
          gdat_data.gdat_params.equil = params_equil.equil;
          gdat_data.equil = equil;
          inb_chord_cxrs=size(gdat_data.data,1);
          inb_time_cxrs=size(gdat_data.data,2);
          psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
          rhopolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
          rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
          rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
          % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
          time_equil=[min(gdat_data.t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.t(end)+0.1)];
          iok=find(~isnan(gdat_data.r(:,1)));
          for itequil=1:length(time_equil)-1
            rr=equil.Rmesh(:,itequil);
            zz=equil.Zmesh(:,itequil);
            psirz_in = equil.psi2D(:,:,itequil);
            it_cxrs_inequil = find(gdat_data.t>time_equil(itequil) & gdat_data.t<=time_equil(itequil+1));
            if ~isempty(it_cxrs_inequil)
              if ~strcmp(upper(gdat_data.gdat_params.source),'CMZ')
                rout=gdat_data.r(iok);
                zout=gdat_data.z(iok);
              else
                rout=gdat_data.r(iok,it_cxrs_inequil);
                zout=gdat_data.z(iok,it_cxrs_inequil);
              end
              psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
              if ~strcmp(upper(gdat_data.gdat_params.source),'CMZ')
                psi_out(iok,it_cxrs_inequil) = repmat(psi_at_routzout,[1,length(it_cxrs_inequil)]);
              else
                psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil));
              end
              rhopolnorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
              for it_cx=1:length(it_cxrs_inequil)
                rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
                rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
              end
            end
          end
          gdat_data.psi = psi_out;
          gdat_data.rhopolnorm = rhopolnorm_out;
          gdat_data.rhotornorm = rhotornorm_out;
          gdat_data.rhovolnorm = rhovolnorm_out;
          aaa.shot = gdat_data.shot;
          aaa.data = gdat_data.data;
          aaa.x = gdat_data.rhopolnorm;
          aaa.t = gdat_data.t;
          aaa.gdat_params = gdat_data.gdat_params;
          aaa = get_grids_1d(aaa,2,1);
          gdat_data.grids_1d = aaa.grids_1d;
          %
          if gdat_data.gdat_params.fit==1
            % add fits
            gdat_data.fit.raw.rhotornorm = NaN*ones(size(gdat_data.ti.data));
            gdat_data.fit.raw.ti.data = NaN*ones(size(gdat_data.ti.data));
            gdat_data.fit.raw.vrot.data = NaN*ones(size(gdat_data.vrot.data));
            rhotornormfit = linspace(0,1,fit_nb_rho_points)';
            gdat_data.fit.rhotornorm = rhotornormfit;
            gdat_data.fit.t = gdat_data.t;
            for it=1:length(gdat_data.t)
              % make rhotor->rhopol transformation for each time since equilibrium might have changed
              gdat_data.fit.rhopolnorm(:,it)=interpos(gdat_data.grids_1d.rhotornorm_equil(:,it), ...
          gdat_data.grids_1d.rhopolnorm_equil(:,it),rhotornormfit);
              idata = find(gdat_data.ti.data(:,it)>0 & gdat_data.rhotornorm(:,it)<1.01);
              if length(idata)>0
                gdat_data.fit.ti.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it);
                gdat_data.fit.ti.raw.data(idata,it) = gdat_data.ti.data(idata,it);
                gdat_data.fit.ti.raw.error_bar(idata,it) = gdat_data.ti.error_bar(idata,it);
                gdat_data.fit.vrot.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it);
                gdat_data.fit.vrot.raw.data(idata,it) = gdat_data.vrot.data(idata,it);
                gdat_data.fit.vrot.raw.error_bar(idata,it) = gdat_data.vrot.error_bar(idata,it);
                [rhoeff,irhoeff] = sort(gdat_data.rhotornorm(idata,it));
                rhoeff = [0; rhoeff];
                % they are some strange low error_bars, so remove these by max(Ti/error_bar)<=10; and changing it to large error bar
                tieff = gdat_data.ti.data(idata(irhoeff),it);
                ti_err_eff = gdat_data.ti.error_bar(idata(irhoeff),it);
                ij=find(tieff./ti_err_eff>10.);
                if ~isempty(ij); ti_err_eff(ij) = tieff(ij)./0.1; end
                vroteff = gdat_data.vrot.data(idata(irhoeff),it);
                vrot_err_eff = gdat_data.vrot.error_bar(idata(irhoeff),it);
                ij=find(vroteff./vrot_err_eff>10.);
                if ~isempty(ij); vrot_err_eff(ij) = vroteff(ij)./0.1; end
                %
                tieff =  [tieff(1); tieff];
                ti_err_eff =  [1e4; ti_err_eff];
                vroteff =  [vroteff(1); vroteff];
                vrot_err_eff =  [1e5; vrot_err_eff];
                [gdat_data.fit.ti.data(:,it), gdat_data.fit.ti.drhotornorm(:,it)] = interpos(rhoeff,tieff,rhotornormfit,fit_tension.ti,[1 0],[0 0],ti_err_eff);
                [gdat_data.fit.vrot.data(:,it), gdat_data.fit.vrot.drhotornorm(:,it)] = interpos(rhoeff,vroteff,rhotornormfit,fit_tension.vrot,[1 0],[0 0],vrot_err_eff);
              end
            end
          end
        end
      end
      if scroll_through_list == 0 || scroll_through_list == 2
        if scroll_through_list == 2
          tmp.(diag_name) = gdat_data;
        end
        if length(gdat_data.gdat_params.source) > i_count
          i_count = i_count + 1;
          diag_name = gdat_data.gdat_params.source{i_count};
        else
          iload = 0;
        end
      elseif scroll_through_list == 1
        if ishotfile_ok == 1
          iload = 0;
        else
          sources_remaining = setdiff(sources_available_for_scroll,sources_tried,'stable');
          if ~isempty(sources_remaining)
            i_count = i_count + 1;
            diag_name = sources_remaining{1};
            sources_tried{i_count} = diag_name;
          else
            iload = 0;
          end
        end
      else
        disp('should not arrive here, check value of scroll_through_list')
        scroll_through_list
        iload = 0;
      end
    end
    if scroll_through_list == 2
      tmp_field=fieldnames(tmp);
      for i=1:length(tmp_field)
        gdat_data.(tmp_field{i}) = tmp.(tmp_field{i});
      end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'ece', 'eced', 'ece_rho', 'eced_rho'}
    nth_points = 13;
    if isfield(gdat_data.gdat_params,'nth_points') && ~isempty(gdat_data.gdat_params.nth_points)
      nth_points = gdat_data.gdat_params.nth_points;
    else
      gdat_data.gdat_params.nth_points = nth_points;
    end
    channels = -1;
    if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels)
      channels = gdat_data.gdat_params.channels;
    end
    if nth_points>=10
      match_rz_to_time = 1;
    else
      match_rz_to_time = 0;
    end
    if isfield(gdat_data.gdat_params,'match_rz_to_time') && ~isempty(gdat_data.gdat_params.match_rz_to_time)
      match_rz_to_time = gdat_data.gdat_params.match_rz_to_time;
    else
      gdat_data.gdat_params.match_rz_to_time = match_rz_to_time;
    end
    time_interval = [];
    if isfield(gdat_data.gdat_params,'time_interval') && ~isempty(gdat_data.gdat_params.time_interval)
      time_interval = gdat_data.gdat_params.time_interval;
    else
      gdat_data.gdat_params.time_interval = time_interval;
    end
    %
    if isfield(gdat_data.gdat_params,'diag_name') && ~isempty(gdat_data.gdat_params.diag_name)
      diag_name = gdat_data.gdat_params.diag_name;
      if strcmp(upper(diag_name),'RMD'); gdat_data.gdat_params.exp_name = 'ECED'; end
    else
      diag_name = 'CEC';
      gdat_data.gdat_params.diag_name = diag_name;
    end
    exp_name_eff = gdat_data.gdat_params.exp_name;
    if strcmp(data_request_eff,'eced')
      exp_name_eff = 'ECED';
      gdat_data.gdat_params.exp_name = exp_name_eff;
    end
    [a,e]=rdaAUG_eff(shot,diag_name,'Trad-A',exp_name_eff,time_interval,gdat_data.gdat_params.extra_arg_sf2sig);
    % [a,e]=rdaAUG_eff(shot,diag_name,'Trad-A',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    inb_chord = size(a.data,1);
    if channels(1)<=0
      channels = [1:inb_chord];
    end
    gdat_data.dim{1} = channels;
    gdat_data.gdat_params.channels = channels;
    if nth_points>1
      gdat_data.data = a.data(channels,1:nth_points:end);
      gdat_data.dim{2} = a.t(1:nth_points:end);
    else
      gdat_data.data = a.data(channels,:);
      gdat_data.dim{2} = a.t;
    end
    gdat_data.x = gdat_data.dim{1};
    gdat_data.t = gdat_data.dim{2};
    gdat_data.dimunits=[{'channels'} ; {'time [s]'}];
    if any(strcmp(fieldnames(a),'units')); gdat_data.units=a.units; end
    try
      [aR,e]=rdaAUG_eff(shot,diag_name,'R-A',exp_name_eff,time_interval,gdat_data.gdat_params.extra_arg_sf2sig);
    catch
    end
    try
      [aZ,e]=rdaAUG_eff(shot,diag_name,'z-A',exp_name_eff,time_interval,gdat_data.gdat_params.extra_arg_sf2sig);
    catch
      disp(['problem with getting z-A in ' diag_name])
    end
    if match_rz_to_time
      % interpolate R structure on ece data time array, to ease plot vs R
      for i=1:length(channels)
        radius.data(i,:) = interp1(aR.t,aR.data(channels(i),:),gdat_data.t);
        zheight.data(i,:) = interp1(aZ.t,aZ.data(channels(i),:),gdat_data.t);
      end
      radiuszheight.t=gdat_data.t;
    else
      radius.data = aR.data(channels,:);
      radiuszheight.t=aR.t;
      zheight.data = aZ.data(channels,:);
    end
    ij=find(gdat_data.data<=0);
    gdat_data.data(ij)=NaN;
    gdat_data.rz.r=radius.data;
    gdat_data.rz.z=zheight.data;
    gdat_data.rz.t = radiuszheight.t;
    gdat_data.data_fullpath = [exp_name_eff '/' diag_name '/Trad-A with R, Z in .r,.z'];

    if strcmp(data_request_eff,'ece_rho') || strcmp(data_request_eff,'eced_rho')
      params_equil = gdat_data.gdat_params;
      params_equil.data_request = 'equil';
      [equil,params_equil,error_status] = gdat_aug(shot,params_equil);
      if error_status>0
        if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end
        return
      end
      gdat_data.gdat_params.equil = params_equil.equil;
      gdat_data.equil = equil;
      gdat_data.data_fullpath = [exp_name_eff '/' diag_name '/Trad-A with .r,.z projected on equil ' gdat_data.gdat_params.equil ' in .rhos'];
      inb_chord_ece=size(gdat_data.rz.r,1);
      inb_time_ece=size(gdat_data.rz.r,2);
      psi_out = NaN*ones(inb_chord_ece,inb_time_ece);
      rhopolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
      rhotornorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
      rhovolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
      % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
      time_equil=[min(gdat_data.rz.t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.rz.t(end)+0.1)];
      for itequil=1:length(time_equil)-1
        rr=equil.Rmesh(:,itequil);
        zz=equil.Zmesh(:,itequil);
        psirz_in = equil.psi2D(:,:,itequil);
        it_ece_inequil = find(gdat_data.rz.t>time_equil(itequil) & gdat_data.rz.t<=time_equil(itequil+1));
        if ~isempty(it_ece_inequil)
          r_it_ece_inequil = gdat_data.rz.r(:,it_ece_inequil);
          iok = find(~isnan(r_it_ece_inequil) & r_it_ece_inequil>0);
          if ~isempty(iok)
            rout=r_it_ece_inequil(iok);
            z_it_ece_inequil = gdat_data.rz.z(:,it_ece_inequil);
            zout=z_it_ece_inequil(iok);
            psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
            [ieff,jeff]=ind2sub(size(gdat_data.rz.r(:,it_ece_inequil)),iok);
            for ij=1:length(iok)
              psi_out(ieff(ij),it_ece_inequil(jeff(ij))) = psi_at_routzout(ij);
            end
            rhopolnorm_out(:,it_ece_inequil) = sqrt(abs((psi_out(:,it_ece_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))));
            for it_cx=1:length(it_ece_inequil)
              rhotornorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
              rhovolnorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
            end
          end
        end
      end
      gdat_data.rhos.psi = psi_out;
      gdat_data.rhos.rhopolnorm = rhopolnorm_out;
      gdat_data.rhos.rhotornorm = rhotornorm_out;
      gdat_data.rhos.rhovolnorm = rhovolnorm_out;
      gdat_data.rhos.t = gdat_data.rz.t;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'eqdsk'}
    %
    time_eqdsks=1.; % default time
    if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time)
      time_eqdsks = gdat_data.gdat_params.time;
    else
      gdat_data.gdat_params.time = time_eqdsks;
      disp(['"time" is expected as an option, choose default time = ' num2str(time_eqdsks)]);
    end
    gdat_data.gdat_params.time = time_eqdsks;
    gdat_data.t = time_eqdsks;
    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;
    dowrite = 1;
    if isfield(gdat_data.gdat_params,'write') && ~isempty(gdat_data.gdat_params.write)
      dowrite = gdat_data.gdat_params.write;
    else
      gdat_data.gdat_params.write = dowrite;
    end
    gdat_data.gdat_params.write = dowrite;
    params_equil = gdat_data.gdat_params;
    params_equil.data_request = 'equil';
    [equil,params_equil,error_status] = gdat_aug(shot,params_equil);
    if error_status>0
      if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end
      return
    end
    gdat_data.gdat_params = params_equil;
    gdat_data.equil = equil;
    if gdat_data.gdat_params.doplot==0
      fignb_handle = -1;
    else
      figure(9999);
      fignb_handle = gcf;
    end
    for itime=1:length(time_eqdsks)
      time_eff = time_eqdsks(itime);
      % use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2
      [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(equil,time_eff,gdat_data.gdat_params.zshift,'source',gdat_data.gdat_params.equil,'extra_arg_sf2sig',gdat_data.gdat_params.extra_arg_sf2sig,'fignb',fignb_handle);
      eqdskAUG.fnamefull = fullfile(['/tmp/' getenv('USER')],['EQDSK_' num2str(shot) 't' num2str(time_eff,'%.4f')]);
      cocos_in = eqdskAUG.cocos;
      if cocos_in ~= equil.cocos;
        equil.cocos = eqdskAUG.cocos;
      end
      if ~isfield(gdat_data.gdat_params,'cocos') || isempty(gdat_data.gdat_params.cocos)
        gdat_data.gdat_params.cocos = cocos_in;
      end
      cocos_out = gdat_data.gdat_params.cocos;
      if cocos_in ~= cocos_out
        [eqdsk_cocosout, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdskAUG,[cocos_in cocos_out]);
      else
        eqdsk_cocosout = eqdskAUG;
      end
      % 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_eqdsks) > 1
        gdat_data.eqdsk{itime} = eqdsk_cocosout;
        if gdat_data.gdat_params.write; gdat_data.eqdsk{itime} = write_eqdsk(eqdskAUG.fnamefull,eqdsk_cocosout); end
        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 = eqdsk_cocosout;
        if gdat_data.gdat_params.write; gdat_data.eqdsk = write_eqdsk(eqdskAUG.fnamefull,eqdsk_cocosout); end
        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 geteqdskAUG with ' gdat_data.gdat_params.equil ' ;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 {'equil'}
    % get equil params and time array in gdat_data.t
    [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data);
    if isempty(Lpf1_t)
      disp('Lpf1_t is empty, probably no data, return')
      return
    end
    % since Lpf depends on time, need to load all first and then loop over time for easier mapping
    [qpsi,e]=rdaAUG_eff(shot,DIAG,'Qpsi',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    ndimrho = size(qpsi.data,2);
    if ndimrho==NTIME_Lpf && ndimrho ~= NTIME
      % data seems to be transposed
      ndimrho = size(qpsi.data,1);
      itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
    else
      itotransposeback = 0;
    end
    qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback);
    ijnan=find(isnan(qpsi.value));
    qpsi.value(ijnan)=0;
    [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
    [phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
    [Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
    [Area,e]=rdaAUG_eff(shot,DIAG,'Area',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    Area=adapt_rda(Area,NTIME,2*ndimrho,itotransposeback);
    [Ri,e]=rdaAUG_eff(shot,DIAG,'Ri',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    Ri=adapt_rda(Ri,NTIME,M_Rmesh,itotransposeback);
    [Zj,e]=rdaAUG_eff(shot,DIAG,'Zj',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    Zj=adapt_rda(Zj,NTIME,N_Zmesh,itotransposeback);
    if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
      [PFM_tree,e]=rdaAUG_eff(shot,DIAG,'PFM',exp_name_eff,[],['''-raw'',' gdat_data.gdat_params.extra_arg_sf2sig]); % -raw necessary for IDE
    else
      [PFM_tree,e]=rdaAUG_eff(shot,DIAG,'PFM',exp_name_eff,[],['''-raw''']); % -raw necessary for IDE
    end
    PFM_tree=adaptPFM_rda(PFM_tree,M_Rmesh,N_Zmesh,NTIME);
    [Pres,e]=rdaAUG_eff(shot,DIAG,'Pres',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    Pres=adapt_rda(Pres,NTIME,2*ndimrho,itotransposeback);
    [Jpol,e]=rdaAUG_eff(shot,DIAG,'Jpol',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    Jpol=adapt_rda(Jpol,NTIME,2*ndimrho,itotransposeback);
    [FFP,e]=rdaAUG_eff(shot,DIAG,'FFP',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    if ~isempty(FFP.value)
      FFP=adapt_rda(FFP,NTIME,ndimrho,itotransposeback);
    else
      FFP.value=NaN*ones(NTIME,max(Lpf1_t));
    end
    if strcmp(DIAG,'EQI') || strcmp(DIAG,'EQH') || strcmp(DIAG,'IDE')
      [Rinv,e]=rdaAUG_eff(shot,DIAG,'Rinv',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      Rinv=adapt_rda(Rinv,NTIME,ndimrho,itotransposeback);
      [R2inv,e]=rdaAUG_eff(shot,DIAG,'R2inv',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      R2inv=adapt_rda(R2inv,NTIME,ndimrho,itotransposeback);
      [Bave,e]=rdaAUG_eff(shot,DIAG,'Bave',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      Bave=adapt_rda(Bave,NTIME,ndimrho,itotransposeback);
      [B2ave,e]=rdaAUG_eff(shot,DIAG,'B2ave',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      B2ave=adapt_rda(B2ave,NTIME,ndimrho,itotransposeback);
      if strcmp(DIAG,'IDE')
        FTRA.value=[];
      else
        [FTRA,e]=rdaAUG_eff(shot,DIAG,'FTRA',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
        FTRA=adapt_rda(FTRA,NTIME,ndimrho,itotransposeback);
      end
    else
      Rinv.value=[]; R2inv.value=[]; Bave.value=[]; B2ave.value=[]; FTRA.value=[];
    end
    [LPFx,e]=rdaAUG_eff(shot,DIAG,'LPFx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    LPFx.value=LPFx.value(1:NTIME); LPFx.data=LPFx.value(1:NTIME); LPFx.t=LPFx.t(1:NTIME);
    [PFxx,e]=rdaAUG_eff(shot,DIAG,'PFxx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    PFxx=adapt_rda(PFxx,NTIME,max(LPFx.value)+1,itotransposeback);
    [RPFx,e]=rdaAUG_eff(shot,DIAG,'RPFx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    RPFx=adapt_rda(RPFx,NTIME,max(LPFx.value)+1,itotransposeback);
    [zPFx,e]=rdaAUG_eff(shot,DIAG,'zPFx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    zPFx=adapt_rda(zPFx,NTIME,max(LPFx.value)+1,itotransposeback);
    % seems "LCFS" q-value is far too large, limit to some max (when diverted)
    max_qValue = 30.; % Note could just put a NaN on LCFS value since ill-defined when diverted
    for it=1:NTIME
      Lpf1 = Lpf1_t(it);
      % Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
      % change it to (radial,time) and use only Lpf+1 points up to LCFS
      ijok=find(abs(qpsi.value)>0); % note: eqr fills in only odd points radially
      % set NaNs to zeroes
      if qpsi.value(ijok(1))<0
        gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
      else
        gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
      end
      % get x values
      gdat_data.psi(:,it)=psi_tree.value(it,Lpf1:-1:1)';
      gdat_data.psi_axis(it)= gdat_data.psi(1,it);
      gdat_data.psi_lcfs(it)= gdat_data.psi(end,it);
      gdat_data.rhopolnorm(:,it) = sqrt(abs((gdat_data.psi(:,it)-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it))));
      if strcmp(DIAG,'EQR');
        % q value has only a few values and from center to edge, assume they are from central rhopol values on
        % But they are every other point starting from 3rd
        ijk=find(gdat_data.qvalue(:,it)~=0);
        if length(ijk)>2
          % now shots have non-zero axis values in eqr
          rhoeff=gdat_data.rhopolnorm(ijk,it);
          qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above
          if ijk(1)>1
            rhoeff = [0.; rhoeff];
            qeff = [qeff(1) ;qeff];
          end
          ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it)));
          qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
          qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]);
        else
          qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
        end
        gdat_data.qvalue(:,it) = qfit;
      end
      % get rhotor values
      gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)';
      gdat_data.rhotornorm(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.phi(end,it)));
      % get rhovol values
      gdat_data.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)';
      gdat_data.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
      gdat_data.rhovolnorm(:,it) = sqrt(abs(gdat_data.vol(:,it) ./ gdat_data.vol(end,it)));
      gdat_data.area(:,it)=Area.value(it,2*Lpf1-1:-2:1)';
      gdat_data.dareadpsi(:,it)=Area.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
      gdat_data.Rmesh(:,it) = Ri.value(it,1:M_Rmesh);
      gdat_data.Zmesh(:,it) = Zj.value(it,1:N_Zmesh);
      gdat_data.psi2D(1:M_Rmesh,1:N_Zmesh,it) = PFM_tree.value(1:M_Rmesh,1:N_Zmesh,it);
      gdat_data.pressure(:,it)=Pres.value(it,2*Lpf1-1:-2:1)';
      gdat_data.dpressuredpsi(:,it)=Pres.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
      if ~isempty(Jpol.value)
        gdat_data.jpol(:,it)=Jpol.value(it,2*Lpf1-1:-2:1)';
        gdat_data.djpolpsi(:,it)=Jpol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
      else
        gdat_data.jpol = [];
        gdat_data.djpolpsi = [];
      end
      if ~strcmp(DIAG,'TRE')
        gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)';
      else
        gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)'/2/pi;
      end
      gdat_data.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1);
      gdat_data.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1);
      gdat_data.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1);
      if ~isempty(Rinv.value)
        gdat_data.rinv(:,it) = Rinv.value(it,Lpf1:-1:1)';
      else
        gdat_data.rinv = [];
      end
      if ~isempty(R2inv.value)
        gdat_data.r2inv(:,it) = R2inv.value(it,Lpf1:-1:1)';
      else
        gdat_data.r2inv = [];
      end
      if ~isempty(Bave.value)
        gdat_data.bave(:,it) = Bave.value(it,Lpf1:-1:1)';
      else
        gdat_data.bave = [];
      end
      if ~isempty(B2ave.value)
        gdat_data.b2ave(:,it) = B2ave.value(it,Lpf1:-1:1)';
      else
        gdat_data.b2ave = [];
      end
      if ~isempty(FTRA.value)
        gdat_data.ftra(:,it) = FTRA.value(it,Lpf1:-1:1)';
      else
        gdat_data.ftra = [];
      end
      %
    end
    gdat_data.x = gdat_data.rhopolnorm;
    %
    try
      [equil_Rcoil,e]=rdaAUG_eff(shot,DIAG,'Rcl',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      gdat_data.Rcoils = equil_Rcoil.value;
      [equil_Zcoil,e]=rdaAUG_eff(shot,DIAG,'Zcl',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      gdat_data.Zcoils = equil_Zcoil.value;
    catch
      gdat_data.Rcoils = [];
      gdat_data.Zcoils = [];
    end
    %
    if strcmp(DIAG,'IDE')
      [IpiPSI,e]=rdaAUG_eff(shot,'IDG','Itor',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      if length(IpiPSI.value)~=NTIME
        disp('Problem with nb time points of IDG/Itor with respect to IDE NTIME, check with O. Sauter')
        return
      end
    else
      [IpiPSI,e]=rdaAUG_eff(shot,DIAG,'IpiPSI',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    end
    gdat_data.Ip = IpiPSI.value(1:NTIME);
    %
    gdat_data.data = gdat_data.qvalue; % put q in data
    gdat_data.units=[]; % not applicable
    gdat_data.data_fullpath = [DIAG ' from exp_name: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
    gdat_data.cocos = 17; % should check FPP (IDE seems 13)
    gdat_data.dim{1} = gdat_data.x;
    gdat_data.dim{2} = gdat_data.t;
    gdat_data.dimunits{1} = 'rhopolnorm';
    gdat_data.dimunits{2} = 'time [s]';

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'gas', 'gas_valve', 'gas_valves'}
    %
    uvs_signals_source = [{'D_tot','Sum of calibrated D2 gas-fluxes (= D/s)'}; ...
                    {'N_tot','Sum of calibrated N2 gas-fluxes (= N*7/s)'}; ...
                    {'H_tot','Sum of calibrated H gas-fluxes (= H/s)'}; ...
                    {'He_tot','Sum of calibrated He gas-fluxes (= He*2/s)'}; ...
                    {'Ne_tot','Sum of calibrated Ne gas-fluxes (= Ne*10/s)'}; ...
                    {'Ar_tot','Sum of calibrated Ar gas-fluxes (= Ar*18/s)'}];
    for i=1:size(uvs_signals_source,1)
      uvs_signals.(uvs_signals_source{i,1}) = gdat_aug(shot,{'UVS',uvs_signals_source{i,1}});
      if isempty(gdat_data.t) && numel(uvs_signals.(uvs_signals_source{i,1}).t)>10
        gdat_data.t = uvs_signals.(uvs_signals_source{i,1}).t; % all have same time base, could use time signal
      end
      if ~isempty(uvs_signals.(uvs_signals_source{i,1}).data) && numel(uvs_signals.(uvs_signals_source{i,1}).data) == numel(gdat_data.t);
        gdat_data.data(i,:) = uvs_signals.(uvs_signals_source{i,1}).data; % all have same time base
      else
        gdat_data.data(i,:) = NaN;
      end
      gdat_data.units{i} = uvs_signals_source{i,2};
    end
    gdat_data.data_fullpath = 'from UVS and uvs_signals_source';
    gdat_data.uvs_signals_source = uvs_signals_source;
    gdat_data.label = uvs_signals_source(:,1)';
    gdat_data.dim{2} = gdat_data.t;
    gdat_data.dimunits{2} = 's';
    gdat_data.dim{1} = [1:size(gdat_data.data,1)];
    gdat_data.dimunits{1} = 'uvs_signals_source index';
    gdat_data.x = gdat_data.dim{1};
    gdat_data.mapping_for.aug.timedim = 1;

    % note also available in UVS, not commented the same way on ISIS
    uvd_signals_source=[{'CFCu01T','PV01','empty comment'}; ...
                        {'CFCo02A','PV02','valve flux, normally for prefill'}; ...
                        {'CFFo02A','PV03','valve flux,'}; ...
                        {'CFA13B','PV04','valve flux, in-vessel pipe for gas imaging'}; ...
                        {'CFA03B','PV05','valve flux,'}; ...
                        {'CFFo07A','PV06','valve flux,'}; ...
                        {'CFDu05X','PV07','valve flux, normally with N2'}; ...
                        {'CFDu05B','PV08','valve flux, normally with D2'}; ...
                        {'CFCo10A','PV09','valve flux, in-vessel pipe into ICRH limiter'}; ...
                        {'CFDu01B','PV10','valve flux, normally with D2'}; ...
                        {'CFDu09B','PV11','valve flux, normally with D2'}; ...
                        {'CFFo10A','PV12','valve flux,'}; ...
                        {'CFFo14A','PV13','valve flux,'}; ...
                        {'CFA03C','PV14','valve flux, short feeding pipe for rare gases'}; ...
                        {'CFA13A','PV15','valve flux,'}; ...
                        {'CFA03A','PV16','valve flux,'}; ...
                        {'CFDu09X','PV17','valve flux, normally with N2'}; ...
                        {'CFDu13X','PV18','valve flux, normally with N2'}; ...
                        {'CFDu13B','PV19','valve flux, normally with D2'}; ...
                        {'CFDu01X','PV20','valve flux, normally with N2'}];
    gdat_data.pvxx.t = [];
    gdat_data.pvxx.t = [];

    for i=1:size(uvd_signals_source,1)
      uvd_signals.(uvd_signals_source{i,2}) = gdat_aug(shot,{'UVD',uvd_signals_source{i,1}});
      if isempty(gdat_data.pvxx.t) && numel(uvd_signals.(uvd_signals_source{i,2}).t)>10
        gdat_data.pvxx.t = uvd_signals.(uvd_signals_source{i,2}).t; % all have same time base
      end
      if ~isempty(uvd_signals.(uvd_signals_source{i,2}).data) && numel(uvd_signals.(uvd_signals_source{i,2}).data) == numel(gdat_data.pvxx.t);
        gdat_data.pvxx.data(i,:) = uvd_signals.(uvd_signals_source{i,2}).data; % all have same time base
      else
        gdat_data.pvxx.data(i,:) = NaN;
      end
      gdat_data.pvxx.units{i} = [uvd_signals_source{i,2} ': ' uvd_signals_source{i,3}];
    end
    gdat_data.pvxx.data_fullpath = 'from UVD and uvd_signals_source';
    gdat_data.pvxx.uvd_signals_source = uvd_signals_source;
    gdat_data.pvxx.label = uvd_signals_source(:,2)';
    gdat_data.pvxx.dim{2} = gdat_data.t;
    gdat_data.pvxx.dimunits{2} = 's';
    gdat_data.pvxx.dim{1} = [1:size(gdat_data.pvxx.data,1)];
    gdat_data.pvxx.dimunits{1} = 'PVxx index';
    gdat_data.pvxx.x = gdat_data.dim{1};

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'ids'}
    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 = [];
      disp('need an ids name in ''source'' parameter');
      disp('check substructure gdat_params.source_available for an ids list');
    end
    ids_gen_ok = exist('ids_gen');
    equil_empty = struct([]);
    if ids_gen_ok ~= 2
      ids_struct_saved = 'ids_structures_20190312.mat';
      if ~exist(ids_struct_saved,'file')
        warning(['function ids_gen not available neither file ids_structures_20190312.mat thus cannot create empty ids: ids_gen(''' ids_top_name ''')']);
        return
      else
        eval(['load ' ids_struct_saved ])
        if isfield(ids_structures,ids_top_name)
          equil_empty = ids_structures.(ids_top_name);
        else
          if ~isempty(ids_top_name);
            warning(['ids ''' ids_top_name ''' does not exist in ids_structures saved in ' ids_struct_saved]);
          end
          gdat_data.gdat_params.source_available = ids_list;
          return
        end
      end
    else
      equil_empty = ids_gen(ids_top_name);
    end
    try
      if ~isempty(shot)
        eval(['[ids_top,ids_top_description]=aug_get_ids_' ids_top_name '(shot,equil_empty);'])
        gdat_data.(ids_top_name) = ids_top;
        gdat_data.([ids_top_name '_description']) = ids_top_description;
      else
        gdat_data.(ids_top_name) = equil_empty;
        gdat_data.([ids_top_name '_description']) = ['shot empty so return default empty structure for ' ids_top_name];
      end
    catch ME_aug_get_ids
      getReport(ME_aug_get_ids)
      disp(['there is a problem with: aug_get_ids_' ids_top_name ...
            ' , may be check if it exists in your path or test it by itself'])
      gdat_data.(ids_top_name) = equil_empty;
      gdat_data.([ids_top_name '_description']) = getReport(ME_aug_get_ids);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'ne','te'}
    exp_name_eff = gdat_data.gdat_params.exp_name;
    % ne or Te from Thomson data on raw z mesh vs (z,t)
    nodenameeff = [upper(data_request_eff(1)) 'e_c'];
    node_child_nameeff = [upper(data_request_eff(1)) 'e_core'];
    [a,error_status]=rdaAUG_eff(shot,'VTA',nodenameeff,exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    if isempty(a.data) || isempty(a.t) || error_status>0
      if gdat_params.nverbose>=3;
        a
        disp(['with data_request= ' data_request_eff])
      end
      return
    end
    gdat_data.(lower(node_child_nameeff)).data = a.data;
    gdat_data.(lower(node_child_nameeff)).t = a.t;
    gdat_data.t = a.t;
    if any(strcmp(fieldnames(a),'units'))
      gdat_data.units=a.units;
    end
    [alow,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'elow_c'],exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    [aup,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'eupp_c'],exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    gdat_data.(lower(node_child_nameeff)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff)).data,gdat_data.(lower(node_child_nameeff)).data-alow.value);
    [a,error_status]=rdaAUG_eff(shot,'VTA','R_core',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    gdat_data.(lower(node_child_nameeff)).r = repmat(a.data,size(gdat_data.(lower(node_child_nameeff)).data,1),1);
    [a,error_status]=rdaAUG_eff(shot,'VTA','Z_core',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    gdat_data.(lower(node_child_nameeff)).z = repmat(a.data',1,size(gdat_data.(lower(node_child_nameeff)).data,2));
    gdat_data.dimunits=[{'Z [m]'} ; {'time [s]'}];
    gdat_data.data_fullpath=[data_request_eff ' from VTA/' upper(data_request_eff(1)) 'e_c and ' upper(data_request_eff(1)) 'e_e'];
    nb_core = size(gdat_data.(lower(node_child_nameeff)).z,1);
    gdat_data.data = [];
    gdat_data.data(1:nb_core,:) = gdat_data.(lower(node_child_nameeff)).data(1:nb_core,:);
    gdat_data.dim=[{gdat_data.(lower(node_child_nameeff)).z};{gdat_data.(lower(node_child_nameeff)).t}];
    gdat_data.error_bar(1:nb_core,:) = gdat_data.(lower(node_child_nameeff)).error_bar(1:nb_core,:);

    % add edge part
    nodenameeff_e = [upper(data_request_eff(1)) 'e_e'];
    node_child_nameeff_e = [upper(data_request_eff(1)) 'e_edge'];
    [a,error_status]=rdaAUG_eff(shot,'VTA',nodenameeff_e,exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    gdat_data.(lower(node_child_nameeff_e)).data = a.data;
    ij_edge_notok = find(a.data>5e21);
    gdat_data.(lower(node_child_nameeff_e)).data(ij_edge_notok) = NaN;
    gdat_data.(lower(node_child_nameeff_e)).t = a.t;
    if ~isempty(a.data)
      [a,error_status]=rdaAUG_eff(shot,'VTA','R_edge',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      gdat_data.(lower(node_child_nameeff_e)).r = repmat(a.data,size(gdat_data.(lower(node_child_nameeff_e)).data,1),1);
      [a,error_status]=rdaAUG_eff(shot,'VTA','Z_edge',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      gdat_data.(lower(node_child_nameeff_e)).z = repmat(a.data',1,size(gdat_data.(lower(node_child_nameeff_e)).data,2));
      nb_edge = size(gdat_data.(lower(node_child_nameeff_e)).z,1);
      iaaa=iround_os(gdat_data.(lower(node_child_nameeff_e)).t,gdat_data.(lower(node_child_nameeff)).t);
      gdat_data.data(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).data(1:nb_edge,iaaa);
      gdat_data.dim{1}(nb_core+1:nb_core+nb_edge,:)=gdat_data.(lower(node_child_nameeff_e)).z(1:nb_edge,iaaa);
      [alow,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'elow_e'],exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      [aup,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'eupp_e'],exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      gdat_data.(lower(node_child_nameeff_e)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff_e)).data, ...
          gdat_data.(lower(node_child_nameeff_e)).data-alow.value);
      gdat_data.error_bar(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).error_bar(1:nb_edge,iaaa);
    else
      gdat_data.(lower(node_child_nameeff_e)).data = [];
      gdat_data.(lower(node_child_nameeff_e)).t = [];
      gdat_data.(lower(node_child_nameeff_e)).error_bar = [];
      gdat_data.(lower(node_child_nameeff_e)).r = [];
      gdat_data.(lower(node_child_nameeff_e)).z = [];
    end
    gdat_data.x=gdat_data.dim{1};

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'ne_rho', 'te_rho', 'nete_rho'}
    sources_available = {'VTA','IDA','TRA'}; % TRA to be included
    if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
      gdat_data.gdat_params.source = 'VTA'; % default source
    end
    sources_exp_name_available = 'AUGD';
    if isfield(gdat_data.gdat_params,'source_exp_name') && ~isempty(gdat_data.gdat_params.source_exp_name)
      sources_exp_name_available = gdat_data.gdat_params.source_exp_name;
    end
    gdat_data.gdat_params.source_exp_name = sources_exp_name_available;
    gdat_data.gdat_params.source = upper(gdat_data.gdat_params.source);
    if strcmp(gdat_data.gdat_params.source(1:2),'EQ');
      warning('set equilibrium choice in ''equil'' parameter and not source, moved to equil')
      gdat_data.gdat_params.equil = gdat_data.gdat_params.source;
      gdat_data.gdat_params.source = [];
    end
    if ~isempty(intersect({'THOMSON'},gdat_data.gdat_params.source)); gdat_data.gdat_params.source = 'VTA'; end
    if ~isempty(intersect({'TRANSP'},gdat_data.gdat_params.source)); gdat_data.gdat_params.source = 'TRA'; end
    if ~isempty(intersect({'IDE'},gdat_data.gdat_params.source)); gdat_data.gdat_params.source = 'IDA'; end
    if isempty(intersect(sources_available,gdat_data.gdat_params.source))
      error(['available source choices: ' sources_available{:} char(10) ' source: ' gdat_data.gdat_params.source ...
             ' not implemented, ask O. Sauter' char(10)]);
    elseif ~strcmp(gdat_data.gdat_params.source,'VTA')
      disp(['At this stage .te, .ne, .fit refer to VTA projected on given equil choice. ' char(10) ...
            'IDA or TRA profiles are added in .ida or .transp respectively']);
    end
    if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit) || ~isnumeric(gdat_data.gdat_params.fit)
      if strcmp(gdat_data.gdat_params.source,'VTA')
        gdat_data.gdat_params.fit = 1; % default do fit (only with VTA source requested)
      else
        % fit profiles filled in by IDA, TRA, ... profiles
        gdat_data.gdat_params.fit = 0;
      end
    end
    %
    params_eff = gdat_data.gdat_params;
    params_eff.data_request=data_request_eff(1:2); % start with ne if nete_rho
    % get raw data
    [gdat_data,params_kin,error_status]=gdat_aug(shot,params_eff);
    gdat_data.gdat_params.data_request=data_request_eff;
    gdat_data.gdat_request=data_request_eff;
    if error_status>0
      if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
      return
    end
    % add rho coordinates
    params_eff.data_request='equil';
    [equil,params_equil,error_status]=gdat_aug(shot,params_eff);
    if error_status>0
      if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
      return
    end
    gdat_data.gdat_params.equil = params_equil.equil;
    %gdat_data.equil = equil;

    % core
    node_child_nameeff = [upper(data_request_eff(1)) 'e_core'];
    inb_chord_core=size(gdat_data.(lower(node_child_nameeff)).r,1);
    inb_time_core=size(gdat_data.(lower(node_child_nameeff)).r,2);
    psi_out_core = NaN*ones(inb_chord_core,inb_time_core);
    rhopolnorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
    rhotornorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
    rhovolnorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
    % edge
    node_child_nameeff_e = [upper(data_request_eff(1)) 'e_edge'];
    inb_chord_edge=size(gdat_data.(lower(node_child_nameeff_e)).r,1);
    inb_time_edge=size(gdat_data.(lower(node_child_nameeff_e)).r,2);
    psi_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
    rhopolnorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
    rhotornorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
    rhovolnorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
    % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
    time_equil=[min(gdat_data.(lower(node_child_nameeff)).t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.(lower(node_child_nameeff)).t(end)+0.1)];
    for itequil=1:length(time_equil)-1
      rr=equil.Rmesh(:,itequil);
      zz=equil.Zmesh(:,itequil);
      psirz_in = equil.psi2D(:,:,itequil);
      it_core_inequil = find(gdat_data.(lower(node_child_nameeff)).t>time_equil(itequil) & gdat_data.(lower(node_child_nameeff)).t<=time_equil(itequil+1));
      if ~isempty(it_core_inequil)
        rout_core=gdat_data.(lower(node_child_nameeff)).r(:,it_core_inequil);
        zout_core=gdat_data.(lower(node_child_nameeff)).z(:,it_core_inequil);
        psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core);
        psi_out_core(:,it_core_inequil) = reshape(psi_at_routzout,inb_chord_core,length(it_core_inequil));
        rhopolnorm_out_core(:,it_core_inequil) = sqrt((psi_out_core(:,it_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
        for it_cx=1:length(it_core_inequil)
          rhotornorm_out_core(:,it_core_inequil(it_cx)) = ...
              interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
          rhovolnorm_out_core(:,it_core_inequil(it_cx)) = ...
              interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
        end
      end
      % edge
      it_edge_inequil = find(gdat_data.(lower(node_child_nameeff_e)).t>time_equil(itequil) & gdat_data.(lower(node_child_nameeff_e)).t<=time_equil(itequil+1));
      if ~isempty(it_edge_inequil)
        rout_edge=gdat_data.(lower(node_child_nameeff_e)).r(:,it_edge_inequil);
        zout_edge=gdat_data.(lower(node_child_nameeff_e)).z(:,it_edge_inequil);
        psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_edge,zout_edge);
        psi_out_edge(:,it_edge_inequil) = reshape(psi_at_routzout,inb_chord_edge,length(it_edge_inequil));
        rhopolnorm_out_edge(:,it_edge_inequil) = sqrt((psi_out_edge(:,it_edge_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
        for it_cx=1:length(it_edge_inequil)
          rhotornorm_out_edge(:,it_edge_inequil(it_cx)) = ...
              interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
          rhovolnorm_out_edge(:,it_edge_inequil(it_cx)) = ...
              interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
        end
      end
    end
    gdat_data.(lower(node_child_nameeff)).psi = psi_out_core;
    gdat_data.(lower(node_child_nameeff)).rhopolnorm = rhopolnorm_out_core;
    gdat_data.(lower(node_child_nameeff)).rhotornorm = rhotornorm_out_core;
    gdat_data.(lower(node_child_nameeff)).rhovolnorm = rhovolnorm_out_core;
    gdat_data.(lower(node_child_nameeff_e)).psi = psi_out_edge;
    gdat_data.(lower(node_child_nameeff_e)).rhopolnorm = rhopolnorm_out_edge;
    gdat_data.(lower(node_child_nameeff_e)).rhotornorm = rhotornorm_out_edge;
    gdat_data.(lower(node_child_nameeff_e)).rhovolnorm = rhovolnorm_out_edge;
    % put values of rhopolnorm for dim{1} by default, all radial mesh for combined core, edge in grids_1d
    gdat_data.x = gdat_data.(lower(node_child_nameeff)).rhopolnorm;
    iaaa=iround_os(gdat_data.(lower(node_child_nameeff_e)).t,gdat_data.(lower(node_child_nameeff)).t);
    gdat_data.x(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhopolnorm(:,iaaa);
    gdat_data.dim{1} = gdat_data.x;
    gdat_data.dimunits{1} = 'rhopolnorm';
    gdat_data2 = get_grids_1d(gdat_data,2,1,gdat_params.nverbose);
    gdat_data.grids_1d = gdat_data2.grids_1d;
% $$$     gdat_data.grids_1d.rhopolnorm = gdat_data.x;
% $$$     gdat_data.grids_1d.psi = gdat_data.(lower(node_child_nameeff)).psi;
% $$$     gdat_data.grids_1d.psi(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).psi(:,iaaa);
% $$$     gdat_data.grids_1d.rhotornorm = gdat_data.(lower(node_child_nameeff)).rhotornorm;
% $$$     gdat_data.grids_1d.rhotornorm(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhotornorm(:,iaaa);
% $$$     gdat_data.grids_1d.rhovolnorm = gdat_data.(lower(node_child_nameeff)).rhovolnorm;
% $$$     gdat_data.grids_1d.rhovolnorm(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhovolnorm(:,iaaa);

    gdat_data.data_fullpath = [gdat_data.data_fullpath ' projected on equilibrium ' gdat_data.gdat_params.equil];

    % if nete_rho, copy data as .ne, get .te and put pe=e ne Te in data:
    gdat_data.(data_request_eff(1:2)).data = gdat_data.data;
    gdat_data.(data_request_eff(1:2)).error_bar = gdat_data.error_bar;
    gdat_data.(data_request_eff(1:2)).units = gdat_data.units;
    gdat_data.(data_request_eff(1:2)).core = gdat_data.([data_request_eff(1:2) '_core']);
    gdat_data.(data_request_eff(1:2)).edge = gdat_data.([data_request_eff(1:2) '_edge']);
    gdat_data = rmfield(gdat_data,{[data_request_eff(1:2) '_core'],[data_request_eff(1:2) '_edge']});
    if strcmp(data_request_eff(1:4),'nete')
      params_eff.data_request=data_request_eff(3:4);
      [gdat_data_te,params_kin,error_status]=gdat_aug(shot,params_eff);
      gdat_data.te.data = gdat_data_te.data;
      gdat_data.te.error_bar = gdat_data_te.error_bar;
      gdat_data.te.units = gdat_data_te.units;
      gdat_data.te.core = gdat_data_te.te_core;
      gdat_data.te.edge = gdat_data_te.te_edge;
      gdat_data.te.error_bar = gdat_data_te.error_bar;
      gdat_data.te.core.psi = gdat_data.ne.core.psi;
      gdat_data.te.core.rhopolnorm = gdat_data.ne.core.rhopolnorm;
      gdat_data.te.core.rhotornorm = gdat_data.ne.core.rhotornorm;
      gdat_data.te.core.rhovolnorm = gdat_data.ne.core.rhovolnorm;
      gdat_data.te.edge.psi = gdat_data.ne.edge.psi;
      gdat_data.te.edge.rhopolnorm = gdat_data.ne.edge.rhopolnorm;
      gdat_data.te.edge.rhotornorm = gdat_data.ne.edge.rhotornorm;
      gdat_data.te.edge.rhovolnorm = gdat_data.ne.edge.rhovolnorm;
      % put pe in main gdat_data
      gdat_data.data = 1.6022e-19.*gdat_data.ne.data.*gdat_data.te.data;
      gdat_data.error_bar = 1.6022e-19 .* (gdat_data.ne.data .* gdat_data.te.error_bar ...
          + gdat_data.te.data .* gdat_data.ne.error_bar);
      gdat_data.units='N/m^2; 1.6022e-19 ne Te';
      gdat_data.data_fullpath=['pe=1.6e-19*ne*Te in data, .ne, .te ' gdat_data.data_fullpath(3:end)];
      gdat_data.label = 'pe';
    end
    % defaults for fits, so user always gets std structure
    gdat_data.fit.rhotornorm = []; % same for both ne and te
    gdat_data.fit.rhopolnorm = [];
    gdat_data.fit.t = [];
    gdat_data.fit.te.data = [];
    gdat_data.fit.te.drhotornorm = [];
    gdat_data.fit.ne.data = [];
    gdat_data.fit.ne.drhotornorm = [];
    gdat_data.fit.raw.te.data = [];
    gdat_data.fit.raw.te.rhotornorm = [];
    gdat_data.fit.raw.ne.data = [];
    gdat_data.fit.raw.ne.rhotornorm = [];
    fit_tension_default = -0.1;
    if isfield(gdat_data.gdat_params,'fit_tension')
      fit_tension = gdat_data.gdat_params.fit_tension;
    else
      fit_tension = fit_tension_default;
    end
    if ~isstruct(fit_tension)
      fit_tension_eff.te = fit_tension;
      fit_tension_eff.ne = fit_tension;
      fit_tension = fit_tension_eff;
    else
      if ~isfield(fit_tension,'te'); fit_tension.te = fit_tension_default; end
      if ~isfield(fit_tension,'ne '); fit_tension.ne = fit_tension_default; end
    end
    gdat_data.gdat_params.fit_tension = fit_tension;
    if isfield(gdat_data.gdat_params,'fit_nb_rho_points')
      fit_nb_rho_points = gdat_data.gdat_params.fit_nb_rho_points;
    else
      fit_nb_rho_points = 201;
    end
    gdat_data.gdat_params.fit_nb_rho_points = fit_nb_rho_points;
    %
    if gdat_data.gdat_params.fit==1
      % add "manual" fits
      gdat_data.fit.t = gdat_data.t;
      rhotornormfit = linspace(0,1,fit_nb_rho_points)';
      gdat_data.fit.rhotornorm = rhotornormfit;
      gdat_data.fit.rhopolnorm = NaN*ones(length(rhotornormfit),length(gdat_data.fit.t));
      if any(strfind(data_request_eff,'te'))
        gdat_data.fit.raw.te.data = NaN*ones(size(gdat_data.te.data));
        gdat_data.fit.raw.te.error_bar = NaN*ones(size(gdat_data.te.data));
        gdat_data.fit.raw.te.rhotornorm = NaN*ones(size(gdat_data.te.data));
        gdat_data.fit.te.data = gdat_data.fit.rhopolnorm;
        gdat_data.fit.te.drhotornorm = gdat_data.fit.rhopolnorm;
      end
      if any(strfind(data_request_eff,'ne'))
        gdat_data.fit.raw.ne.data = NaN*ones(size(gdat_data.ne.data));
        gdat_data.fit.raw.ne.error_bar = NaN*ones(size(gdat_data.ne.data));
        gdat_data.fit.raw.ne.rhotornorm = NaN*ones(size(gdat_data.ne.data));
        gdat_data.fit.ne.data =gdat_data.fit.rhopolnorm;
        gdat_data.fit.ne.drhotornorm = gdat_data.fit.rhopolnorm;
      end
      for it=1:length(gdat_data.t)
        % make rhotor->rhopol transformation for each time since equilibrium might have changed
        gdat_data.fit.rhopolnorm(:,it)=interpos(gdat_data.grids_1d.rhotornorm_equil(:,it), ...
          gdat_data.grids_1d.rhopolnorm_equil(:,it),rhotornormfit);
        gdat_data.fit.rhopolnorm(:,it)=interpos(gdat_data.grids_1d.rhotornorm_equil(:,it), ...
          gdat_data.grids_1d.rhopolnorm_equil(:,it),rhotornormfit);
        if any(strfind(data_request_eff,'te'))
          idatate = find(gdat_data.te.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05);
          if length(idatate)>0
            gdat_data.fit.raw.te.rhotornorm(idatate,it) = gdat_data.grids_1d.rhotornorm(idatate,it);
            gdat_data.fit.raw.te.data(idatate,it) = gdat_data.te.data(idatate,it);
            gdat_data.fit.raw.te.error_bar(idatate,it) = gdat_data.te.error_bar(idatate,it);
            [rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhotornorm(idatate,it));
            rhoeff = [0; rhoeff];
            teeff = gdat_data.te.data(idatate(irhoeff),it);
            te_err_eff = gdat_data.te.error_bar(idatate(irhoeff),it);
            % they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
            ij=find(teeff./te_err_eff>10.);
            if ~isempty(ij); te_err_eff(ij) = max(max(te_err_eff),teeff(ij)./0.1); end
            %
            teeff =  [teeff(1); teeff];
            te_err_eff =  [1e4; te_err_eff];
            % add some points to reach 0 at 1.05
            if max(rhoeff)<1.03
              rhoeff = [rhoeff; 1.05];
              teeff = [teeff; 0.];
              te_err_eff(end) = 10.*max(te_err_eff(2:end));
              te_err_eff = [te_err_eff; min(te_err_eff)]; % to give more weight for new last point
            else
              teeff(end) = 0.;
              te_err_eff(end-1) = 10.*max(te_err_eff(2:end));
              te_err_eff(end) = min(te_err_eff); % to give more weight for new value of last point
            end
            [gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhotornorm(:,it)] = interpos(rhoeff,teeff,rhotornormfit,fit_tension.te,[1 0],[0 0],te_err_eff);
            % avoid neg points or positive edge gradient
            if ~isempty(find(gdat_data.fit.te.data(:,it)<0)) || gdat_data.fit.te.drhotornorm(end,it) > 0
              ij= find(rhoeff>=0.85 & rhoeff<1.04);
              if ~isempty(ij)
                te_err_eff(ij) = 100.*min(te_err_eff);
              else
                te_err_eff(end-min(5,length(te_err_eff)-3):end-1) = 100.*min(te_err_eff);
              end
              [gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhotornorm(:,it)] = interpos(rhoeff,teeff,rhotornormfit,fit_tension.te,[1 0],[0 0],te_err_eff);
            end
          end
        end
        if any(strfind(data_request_eff,'ne'))
          idatane = find(gdat_data.ne.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05);
          if length(idatane)>0
            gdat_data.fit.raw.ne.rhotornorm(idatane,it) = gdat_data.grids_1d.rhotornorm(idatane,it);
            gdat_data.fit.raw.ne.data(idatane,it) = gdat_data.ne.data(idatane,it);
            gdat_data.fit.raw.ne.error_bar(idatane,it) = gdat_data.ne.error_bar(idatane,it);
            [rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhotornorm(idatane,it));
            rhoeff = [0; rhoeff];
            % they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
            neeff = gdat_data.ne.data(idatane(irhoeff),it);
            ne_err_eff = gdat_data.ne.error_bar(idatane(irhoeff),it);
            ij=find(neeff./ne_err_eff>10.);
            if ~isempty(ij); ne_err_eff(ij) = max(max(ne_err_eff),neeff(ij)./0.1); end
            %
            neeff =  [neeff(1); neeff];
            ne_err_eff =  [1e21; ne_err_eff];
            % add some points to reach 0 at 1.05
            if max(rhoeff)<1.03
              rhoeff = [rhoeff; 1.05];
              neeff = [neeff; 0.];
              ne_err_eff(end) = 10.*max(ne_err_eff(2:end)); % to give more weight for new last point
              ne_err_eff = [ne_err_eff; min(ne_err_eff)];
            else
              neeff(end) = 0.;
              ne_err_eff(end-1) = 10.*max(ne_err_eff(2:end));
              ne_err_eff(end) = min(ne_err_eff); % to give more weight for new value of last point
            end
            [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhotornorm(:,it)] = interpos(rhoeff,neeff,rhotornormfit,fit_tension.ne,[1 0],[0 0],ne_err_eff);
            % avoid neg points or positive edge gradient
            if ~isempty(find(gdat_data.fit.ne.data(:,it)<0)) || gdat_data.fit.ne.drhotornorm(end,it) > 0
              ij= find(rhoeff>=0.85 & rhoeff<1.04);
              if ~isempty(ij)
                ne_err_eff(ij) = 100.*min(ne_err_eff);
              else
                ne_err_eff(end-min(5,length(ne_err_eff)-3):end-1) = 100.*min(ne_err_eff);
              end
              [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhotornorm(:,it)] = interpos(rhoeff,neeff,rhotornormfit,fit_tension.ne,[1 0],[0 0],ne_err_eff);
            end
          end
        end
      end
    end
    % add other profiles in fit (without raw data part)
    switch gdat_data.gdat_params.source
     case 'IDA'
      gdat_data.fit.full_path = 'from IDA: Te, Te_unc, ne, ne_unc, rhop and rhot';
      params_eff.data_request={'IDA','','','area-base:rhop'};
      [gdat_data_ida_rhop,params_kin,error_status]=gdat_aug(shot,params_eff);
      gdat_data.fit.rhopolnorm = gdat_data_ida_rhop.data;
      params_eff.data_request={'IDA','rhot'}; % it is actually rhotor_norm
      [gdat_data_te_ida_rhot,params_kin,error_status]=gdat_aug(shot,params_eff);
      gdat_data.fit.rhotornorm = gdat_data_te_ida_rhot.data;
      if strcmp(data_request_eff(1:4),'nete') || strcmp(data_request_eff(1:2),'te')
        % Te
        params_eff.data_request={'IDA','Te'};
        [gdat_data_te_ida,params_kin,error_status]=gdat_aug(shot,params_eff);
        gdat_data.fit.te.data = gdat_data_te_ida.data;
        gdat_data.fit.t = gdat_data_te_ida.t;
        params_eff.data_request={'IDA','Te_unc'};
        [gdat_data_te_ida_err,params_kin,error_status]=gdat_aug(shot,params_eff);
        gdat_data.fit.te.error_bar = gdat_data_te_ida_err.data;
      end
      if strcmp(data_request_eff(1:4),'nete') || strcmp(data_request_eff(1:2),'ne')
        % ne
        params_eff.data_request={'IDA','ne'};
        [gdat_data_ne_ida,params_kin,error_status]=gdat_aug(shot,params_eff);
        gdat_data.fit.ne.data = gdat_data_ne_ida.data;
        gdat_data.fit.t = gdat_data_ne_ida.t;
        params_eff.data_request={'IDA','ne_unc'};
        [gdat_data_ne_ida_err,params_kin,error_status]=gdat_aug(shot,params_eff);
        gdat_data.fit.ne.error_bar = gdat_data_ne_ida_err.data; % time as 2nd dim
      end
     case 'TRA'
      gdat_data.fit.full_path = 'from TRA: TE, NE, rhop from PLFLX and XB (.x is rhotornorm mid-points, XB end-point)';
      % (tra.XB.data(1:end,it)'+[0 ,tra.XB.data(1:end-1,it)'])/2 = tra.TE.x(1:end,it)'
      if strcmp(data_request_eff(1:4),'nete') || strcmp(data_request_eff(1:2),'te')
        % Te
        params_eff.data_request={'TRA','TE',gdat_data.gdat_params.source_exp_name};
        [gdat_data_te_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
        gdat_data.fit.te.data = gdat_data_te_tra.data;
        gdat_data.fit.t = gdat_data_te_tra.t;
        gdat_data.fit.rhotornorm = gdat_data_te_tra.x;
        params_eff.data_request={'TRA','PLFLX',gdat_data.gdat_params.source_exp_name};
        [gdat_data_plflx_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
        rhopol_endpoints = sqrt((gdat_data_plflx_tra.data-0)./(repmat(gdat_data_plflx_tra.data(end,:),size(gdat_data_plflx_tra.data,1),1)-0.));
        for it=1:size(gdat_data_plflx_tra.data,2)
          gdat_data.fit.rhopolnorm(:,it) = interpos([0. ; gdat_data_plflx_tra.x(:,it)],[0. ; rhopol_endpoints(:,it)],gdat_data.fit.rhotornorm(:,it));
        end
      end
      if strcmp(data_request_eff(1:4),'nete') || strcmp(data_request_eff(1:2),'ne')
        % ne
        params_eff.data_request={'TRA','NE',gdat_data.gdat_params.source_exp_name};
        [gdat_data_ne_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
        gdat_data.fit.ne.data = gdat_data_ne_tra.data .* 1e6; % ne in cm^3 !!
        if ~strcmp(data_request_eff(1:4),'nete')
          gdat_data.fit.t = gdat_data.t;
          gdat_data.fit.rhotornorm = gdat_data_ne_tra.x;
          params_eff.data_request={'TRA','PLFLX',gdat_data.gdat_params.source_exp_name};
          [gdat_data_plflx_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
          rhopol_endpoints = sqrt((gdat_data_plflx_tra.data-0)./(repmat(gdat_data_plflx_tra.data(end,:),size(gdat_data_plflx_tra.data,1),1)-0.));
          for it=1:size(gdat_data_plflx_tra.data,2)
            gdat_data.fit.rhopolnorm(:,it) = interpos([0. ; gdat_data_plflx_tra.x(:,it)],[0. ; rhopol_endpoints(:,it)],gdat_data.fit.rhotornorm(:,it));
          end
        end
      end
     otherwise
      if ~strcmp(gdat_data.gdat_params.source,'VTA')
        disp(['profiles ''fit'' from ' gdat_data.gdat_params.source ' not defined. At this stage only: ' sources_available{:} ...
              ' are available. Ask O. Sauter if need be']);
      end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'pgyro'}
    extra_arg_sf2sig_eff_string = '';
    if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
      extra_arg_sf2sig_eff_string = [',' gdat_data.gdat_params.extra_arg_sf2sig];
    end
    %  LOAD MULTI CHANNEL DATA ECS
    %  powers, frequencies, etc
    params_eff = gdat_data.gdat_params;
    params_eff.data_request={'ECS','PECRH'};
    % pgyro tot in index=9
    try
      gdat_data=gdat_aug(shot,params_eff);
      gdat_data.data_request = data_request_eff;
      gdat_data.gdat_params.data_request = data_request_eff;
    catch
      if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
      gdat_data.data_request = data_request_eff;
      return
    end
    nb_timepoints = length(gdat_data.t);
    pgyro = NaN*ones(nb_timepoints,9);
    pgyro(:,9) = reshape(gdat_data.data,nb_timepoints,1);
    gdat_data.data = pgyro;
    labels{9} = 'ECtot';
    for i=1:4
      % "old" ECRH1 gyrotrons: gyro 1 to 4 in pgyro
      params_eff.data_request={'ECS',['PG' num2str(i)]};
      gdat_data_i=gdat_aug(shot,params_eff);
      if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
      else
        gdat_data.ec{i} = gdat_data_i;
        gdat_data.data(:,i) = reshape(gdat_data_i.data,nb_timepoints,1);
        gdat_data.dim = [{gdat_data_i.t} {[1:9]}];
        if max(gdat_data_i.data) > 0.
          labels{i} = ['EC_' num2str(i)];
        end
      end
      sys1_source = 'P_sy1_g';
      if shot > 33725; sys1_source = 'P_sy3_g'; end
      try
        [a,e]=rdaAUG_eff(shot,'ECS',[sys1_source num2str(i)],gdat_data.gdat_params.exp_name,[], ...
          gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'gyr_freq']);
        % eval(['a = sf2par(''ECS'',shot,''gyr_freq'','sys1_source num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
      catch
        % gyr_freq not present (old shots for example)
        a=[];
      end
      if isempty(a)
      else
        gdat_data.ec{i}.freq = a;
        gdat_data.freq_ec(i) = a.value;
      end
      try
        [a,e]=rdaAUG_eff(shot,'ECS',[sys1_source num2str(i)],gdat_data.gdat_params.exp_name,[], ...
          gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GPolPos']);
        % eval(['a = sf2par(''ECS'',shot,''GPolPos'','sys1_source num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
      catch
        % GPolPos not present
        a=[];
      end
      if isempty(a)
      else
        gdat_data.ec{i}.polpos = a.value;
        gdat_data.polpos_ec(i) = a.value;
      end
      try
        [a,e]=rdaAUG_eff(shot,'ECS',[sys1_source num2str(i)],gdat_data.gdat_params.exp_name,[], ...
          gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GTorPos']);
        % eval(['a = sf2par(''ECS'',shot,''GTorPos'','sys1_source num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
      catch
        a=[];
      end
      if isempty(a)
      else
        gdat_data.ec{i}.torpos = a.value;
        gdat_data.torpos_ec(i) = a.value;
      end
      % "new" ECRH2 gyrotrons: gyro 5 to 8 in pgyro
      params_eff.data_request={'ECS',['PG' num2str(i) 'N']};
      gdat_data_i=gdat_aug(shot,params_eff);
      if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
      else
        gdat_data.ec{i+4} = gdat_data_i;
        gdat_data.data(:,i+4) = reshape(gdat_data_i.data,nb_timepoints,1);
        gdat_data.dim = [{gdat_data_i.t} {[1:9]}];
        if max(gdat_data_i.data) > 0.
          labels{i+4} = ['EC_' num2str(i+4)];
        end
      end
      try
        [a,e]=rdaAUG_eff(shot,'ECS',['P_sy2_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ...
          gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'gyr_freq']);
        % eval(['a = sf2par(''ECS'',shot,''gyr_freq'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
      catch
        a=[];
      end
      if isempty(a)
      else
        gdat_data.ec{i+4}.freq = a;
        gdat_data.freq_ec(i+4) = a.value;
      end
      try
        [a,e]=rdaAUG_eff(shot,'ECS',['P_sy2_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ...
          gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GPolPos']);
        % eval(['a = sf2par(''ECS'',shot,''GPolPos'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
      catch
        a=[];
      end
      if isempty(a)
      else
        gdat_data.ec{i+4}.polpos = a.value;
        gdat_data.polpos_ec(i+4) = a.value;
      end
      try
        [a,e]=rdaAUG_eff(shot,'ECS',['P_sy2_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ...
          gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GTorPos']);
        % eval(['a = sf2par(''ECS'',shot,''GTorPos'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
      catch
        a=[];
      end
      if isempty(a)
      else
        gdat_data.ec{i+4}.torpos = a.value;
        gdat_data.torpos_ec(i+4) = a.value;
      end
      params_eff.data_request={'ECN',['G' num2str(i) 'POL']};
      gdat_data_i=gdat_aug(shot,params_eff);
      if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
      else
        gdat_data.ec{i+4}.gpol_ec = gdat_data_i;
      end
      params_eff.data_request={'ECN',['G' num2str(i) 'TOR']};
      gdat_data_i=gdat_aug(shot,params_eff);
      if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
      else
        gdat_data.ec{i+4}.gtor_ec = gdat_data_i;
      end
      ishift=0;
      if shot >= 36794; ishift=4; end
      params_eff.data_request={'ECN',['G' num2str(i+ishift) 'PO4']};
      try
        gdat_data_i=gdat_aug(shot,params_eff);
      catch ME1
        disp(getReport(ME1))
        gdat_data_i.data = [];
      end
      if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
      else
        gdat_data.ec{i+4}.gpo4_ec = gdat_data_i;
      end
      params_eff.data_request={'ECN',['G' num2str(i+ishift) 'PO8']};
      try
        gdat_data_i=gdat_aug(shot,params_eff);
      catch ME2
        disp(getReport(ME2))
        gdat_data_i.data = [];
      end
      if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
      else
        gdat_data.ec{i+4}.gpo8_ec = gdat_data_i;
      end
    end
    if ~isempty(gdat_data.dim)
      gdat_data.t = gdat_data.dim{1};
      gdat_data.x = gdat_data.dim{2};
      gdat_data.dimunits=[{'time [s]'} {'ECRH1(1:4) ECRH2(1:4) ECtot'}];
      gdat_data.units='W';
      gdat_data.freq_ech_units = 'GHz';
      gdat_data.data_fullpath=['ECS/' 'PGi and PGiN, etc'];
      icount=0;
      for i=1:length(labels)
        if ~isempty(labels{i})
          icount=icount+1;
          gdat_data.label{icount} = labels{i};
        end
      end
    else
      gdat_data.freq_ech_units =[]';
      gdat_data.ec = [];
      gdat_data.freq_ec = [];
      gdat_data.polpos_ec = [];
      gdat_data.torpos_ec = [];
    end


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'powers'}
    sources_avail = {'ohm','ec','nbi','ic','rad'}; % note should allow ech, nbh, ohmic in parameter sources
    for i=1:length(sources_avail)
      gdat_data.(sources_avail{i}).data = [];
      gdat_data.(sources_avail{i}).units = [];
      gdat_data.(sources_avail{i}).dim=[];
      gdat_data.(sources_avail{i}).dimunits=[];
      gdat_data.(sources_avail{i}).t=[];
      gdat_data.(sources_avail{i}).x=[];
      gdat_data.(sources_avail{i}).data_fullpath=[];
      gdat_data.(sources_avail{i}).label=[];
    end
    if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
      gdat_data.gdat_params.source = sources_avail;
    elseif ~iscell(gdat_data.gdat_params.source)
      if ischar(gdat_data.gdat_params.source)
        gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source);
        if ~any(strmatch(gdat_data.gdat_params.source,lower(sources_avail)))
          if (gdat_params.nverbose>=1)
            warning(['source= ' gdat_data.gdat_params.source ' is not part of the available sources: ' sprintf('''%s'' ',sources_avail{:})]);
          end
          return
        else
          gdat_data.gdat_params.source = {gdat_data.gdat_params.source};
        end
      else
        if (gdat_params.nverbose>=1); warning([' source parameter not compatible with: ' sprintf('''%s'' ',sources_avail{:})]); end
        return
      end
    else
      for i=1:length(gdat_data.gdat_params.source)
        gdat_data.gdat_params.source{i} = lower(gdat_data.gdat_params.source{i});
        if ~any(strmatch(gdat_data.gdat_params.source{i},lower(sources_avail)))
          if gdat_data.gdat_params.nverbose>=1
            warning(['source = ' gdat_data.gdat_params.source{i} ' not expected with data_request= ' data_request_eff])
          end
        end
      end
    end
    % always start from ohmic so can use this time as base time since should yield full shot

    fields_to_copy = {'data','units','dim','dimunits','t','x','data_fullpath','label','help','gdat_params'};
    fields_to_not_copy = {'shot','gdat_request'};
    % total of each source in .data, but full data in subfield like pgyro in .ec, to check for nbi
    params_eff = gdat_data.gdat_params;
    % ohmic, use its time-base
    params_eff.data_request={'TOT','P_OH'};
    try
      ohm=gdat_aug(shot,params_eff);
    catch
      ohm.data = [];
      ohm.dim = [];
    end
    if ~isempty(ohm.data) && ~isempty(ohm.dim)
      for i=1:length(fields_to_copy)
        if isfield(ohm,fields_to_copy{i})
          gdat_data.ohm.(fields_to_copy{i}) = ohm.(fields_to_copy{i});
        end
      end
      gdat_data.ohm.raw_data = gdat_data.ohm.data;
    else
      if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
      return
    end
    mapping_for_aug.timedim = 1; mapping_for_aug.gdat_timedim = 1;
    taus = -10;
    %
    % add each source in main.data, on ohm time array
    gdat_data.t = gdat_data.ohm.t;
    gdat_data.dim{1} = gdat_data.t;
    gdat_data.dimunits{1} = 's';
    gdat_data.ohm.data = interpos(gdat_data.t,gdat_data.ohm.raw_data,5.*taus);
    gdat_data.data = reshape(gdat_data.ohm.data,length(gdat_data.t),1);
    gdat_data.ohm.tension = 5.*taus;
    gdat_data.x =[1];
    gdat_data.label={'P_{ohm}'};
    gdat_data.units = 'W';
    %
    if any(strmatch('ec',gdat_data.gdat_params.source))
      % ec
      params_eff.data_request={'ECS','PECRH'};
      params_eff.data_request='pgyro';
      try
        ec=gdat_aug(shot,params_eff);
      catch
      end
      if ~isempty(ec.data) && ~isempty(ec.dim)
        for i=1:length(fields_to_copy)
          % if has pgyro, use not_copy
          if isfield(ec,fields_to_copy{i}) && ~any(strmatch(fields_to_not_copy,fields_to_copy{i}))
            gdat_data.ec.(fields_to_copy{i}) = ec.(fields_to_copy{i});
          end
        end
        gdat_data.data(:,end+1) = interpos(-21,gdat_data.ec.t,gdat_data.ec.data(:,end),gdat_data.t);
        gdat_data.x(end+1) =gdat_data.x(end)+1;
        gdat_data.label{end+1}='P_{ec}';
      end
    end
    %
    if any(strmatch('nb',gdat_data.gdat_params.source))
      % nbi
      params_eff.data_request={'NIS','PNIQ'};
      try
        nbi=gdat_aug(shot,params_eff);
      catch
        nbi.data = [];
        nbi.dim = [];
      end
      if ~isempty(nbi.data) && ~isempty(nbi.dim)
        for i=1:length(fields_to_copy)
          if isfield(nbi,fields_to_copy{i})
            gdat_data.nbi.(fields_to_copy{i}) = nbi.(fields_to_copy{i});
          end
        end
        % add to main with linear interpolation and 0 for extrapolated values
        gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi.t,gdat_data.nbi.data(:,end),gdat_data.t);
        gdat_data.x(end+1) =gdat_data.x(end)+1;
        gdat_data.label{end+1}='P_{nbi}';
      end
    end
    %
    if any(strmatch('ic',gdat_data.gdat_params.source))
      % ic
      params_eff.data_request={'ICP','PICRN'};
      try
        ic=gdat_aug(shot,params_eff);
      catch
        ic.data = [];
        ic.dim = [];
      end
      if ~isempty(ic.data) && ~isempty(ic.dim)
        for i=1:length(fields_to_copy)
          if isfield(ic,fields_to_copy{i})
            gdat_data.ic.(fields_to_copy{i}) = ic.(fields_to_copy{i});
          end
        end
        gdat_data.data(:,end+1) = interpos(-21,gdat_data.ic.t,gdat_data.ic.data,gdat_data.t);
        gdat_data.x(end+1) =gdat_data.x(end)+1;
        gdat_data.label{end+1}='P_{ic}';
      end
    end
    index_prad = [];
    if any(strmatch('rad',gdat_data.gdat_params.source))
      % radiated power
      params_eff.data_request='prad';
      try
        rad=gdat_aug(shot,params_eff);
      catch
      end
      if ~isempty(rad.data) && ~isempty(rad.dim)
        for i=1:length(fields_to_copy)
          if isfield(rad,fields_to_copy{i})
            gdat_data.rad.(fields_to_copy{i}) = rad.(fields_to_copy{i});
          end
        end
        % add to main with linear interpolation and 0 for extrapolated values
        gdat_data.data(:,end+1) = interpos(-21,gdat_data.rad.t,gdat_data.rad.data,gdat_data.t);
        index_prad = size(gdat_data.data,2); % to not add for ptot
        gdat_data.x(end+1) =gdat_data.x(end)+1;
        gdat_data.label{end+1}='P_{rad}';
      end
    end
    % add tot power
    index_noprad = setdiff([1:size(gdat_data.data,2)],index_prad);
    gdat_data.data(:,end+1) = sum(gdat_data.data(:,index_noprad),2);
    gdat_data.label{end+1}='P_{tot,h}';
    gdat_data.x(end+1) =gdat_data.x(end)+1;
    gdat_data.dim{2} = gdat_data.x; gdat_data.dimunits{2} = '';
    gdat_data.data_fullpath = 'tot powers from each sources, and total heating power in .data(:,end)';

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'q_rho'}
    [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data);
    % since Lpf depends on time, need to load all first and then loop over time for easier mapping
    [qpsi,e]=rdaAUG_eff(shot,DIAG,'Qpsi',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    ndimrho = size(qpsi.data,2);
    if ndimrho==NTIME_Lpf
      % data seems to be transposed
      ndimrho = size(qpsi.data,1);
      itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
    else
      itotransposeback = 0;
    end
    qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback);
    ijnan=find(isnan(qpsi.value));
    qpsi.value(ijnan)=0;
    [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
    [phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
    [Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
    % seems "LCFS" q-value is far too large, limit to some max (when diverted)
    max_qValue = 30.; % Note could just put a NaN on LCFS value since ill-defined when diverted
    for it=1:NTIME
      Lpf1 = Lpf1_t(it);
      % Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
      % change it to (radial,time) and use only Lpf+1 points up to LCFS
      ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially
      % set NaNs to zeroes
      if qpsi.value(ijok(1),1)<0
        gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
      else
        gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
      end
      % get x values
      psi_it=psi_tree.value(it,Lpf1:-1:1)';
      gdat_data.psi_axis(it)= psi_it(1);
      gdat_data.psi_lcfs(it)= psi_it(end);
      gdat_data.rhopolnorm(:,it) = sqrt(abs((psi_it-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it))));
      if strcmp(DIAG,'EQR');
        % q value has only a few values and from center to edge, assume they are from central rhopol values on
        % But they are every other point starting from 3rd
        ijk=find(gdat_data.qvalue(:,it)~=0);
        if length(ijk)>2
          % now shots have non-zero axis values in eqr
          rhoeff=gdat_data.rhopolnorm(ijk,it);
          qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above
          if ijk(1)>1
            rhoeff = [0.; rhoeff];
            qeff = [qeff(1) ;qeff];
          end
          ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it)));
          qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
          qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]);
        else
          qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
        end
        gdat_data.qvalue(:,it) = qfit;
      end
      % get rhotor values
      phi_it = phi_tree.value(it,Lpf1:-1:1)';
      gdat_data.rhotornorm(:,it) = sqrt(abs(phi_it ./ phi_it(end)));
      % get rhovol values
      vol_it=Vol.value(it,2*Lpf1-1:-2:1)';
      gdat_data.rhovolnorm(:,it) = sqrt(abs(vol_it ./ vol_it(end)));
      % Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
      % change it to (radial,time) and use only Lpf+1 points up to LCFS
      ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially
      % max value 1e6, change to 2*extrapolation (was max_qValue before)
      q_edge=interpos(gdat_data.rhotornorm(1:end-1,it),qpsi.value(it,Lpf1:-1:2),1,-0.1);
      gdat_data.qvalue(:,it) = qpsi.value(it,Lpf1:-1:1)';
      if abs(gdat_data.qvalue(end,it)) > 1e3
        % assume diverted
        gdat_data.qvalue(end,it) = 2. * q_edge;
      end
% $$$       if qpsi.value(ijok(1),1)<0
% $$$   gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
% $$$       else
% $$$   gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
% $$$       end
    end
    gdat_data.x = gdat_data.rhopolnorm;
    % get time values
    gdat_data.data = gdat_data.qvalue; % put q in data
    gdat_data.units=[]; % not applicable
    gdat_data.data_fullpath = [DIAG ' from exp_name: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
    gdat_data.cocos = 17; % should check FPP (IDE seems 13)
    gdat_data.dim{1} = gdat_data.x;
    gdat_data.dim{2} = gdat_data.t;
    gdat_data.dimunits{1} = 'rhopolnorm';
    gdat_data.dimunits{2} = 'time [s]';
    % add grids_1d to have rhotor, etc
    gdat_data = get_grids_1d(gdat_data,2,1,gdat_params.nverbose);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'psi_axis', 'psi_edge'}
    if strcmp(upper(gdat_data.gdat_params.equil),'FPG')
      gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev;
      gdat_params_eff.data_request = [{'FPG'},{'fax-bnd'},{'AUGD'}];
      gdat_data = gdat_aug(gdat_data.shot,gdat_params_eff);
      gdat_data.gdat_params = gdat_params_prev;
      gdat_data.label = 'psi\_axis-psi\_edge';
      gdat_data.data_fullpath = [gdat_data.data_fullpath ' yields only psi\_axis-psi\_edge from FPG/fax-bnd'];
    else
      [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data);
      % since Lpf depends on time, need to load all first and then loop over time for easier mapping
      [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
      ndimrho = size(psi_tree.data,2);
      if ndimrho==NTIME_Lpf
        % data seems to be transposed
        ndimrho = size(psi_tree.data,1);
        itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
      else
        itotransposeback = 0;
      end
      psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
      ijnan=find(isnan(psi_tree.value));
      psi_tree.value(ijnan)=0;
      gdat_data.dim{1} = gdat_data.t;
      gdat_data.dimunits{1} = 'time [s]';
      for it=1:NTIME
        Lpf1 =min(Lpf1_t(it),size(psi_tree.value,2));
        psi_it=psi_tree.value(it,Lpf1:-1:1)';
        if strcmp(data_request_eff,'psi_axis')
          gdat_data.data(it)= psi_it(1);
        elseif strcmp(data_request_eff,'psi_edge')
          gdat_data.data(it)= psi_it(end);
        else
        end
      end
      gdat_data.units = psi_tree.units;
      gdat_data.data_fullpath = [DIAG '/PFL extract ' data_request_eff];
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'raptor'}
    sources_avail = {'observer','predictive'};
    for i=1:length(sources_avail)
      gdat_data.(sources_avail{i}).data = [];
      gdat_data.(sources_avail{i}).units = [];
      gdat_data.(sources_avail{i}).dim=[];
      gdat_data.(sources_avail{i}).dimunits=[];
      gdat_data.(sources_avail{i}).t=[];
      gdat_data.(sources_avail{i}).x=[];
      gdat_data.(sources_avail{i}).data_fullpath=[];
      gdat_data.(sources_avail{i}).label=[];
    end
    rap_signals={'ECE_f','Hmode','Ip','P_LH','Phib','beta','cs','exflag','g_f','it','li3','ne_f','res','tcomp','ECEm', ...
                 'F','Vp','chie','dk_k','dpsi','dte','g1','g23or','g2','g3','iota','jbs','jec','jnb','jpar','ne','ni', ...
                 'p','paux','pbrem','pec','pei','pnbe','poh','prad','psi','q','rhtECE','sdk_k','shear','signeo','sxk_k', ...
                 'te','ti','upl','xk_k','ze'};

    rap_desc_signals={'rts:Dia/RAPTOR/ECEflag.val','rts:Dia/RAPTOR/Hmode.val','rts:Dia/RAPTOR/Ip.val','rts:Dia/RAPTOR/P_LH.val',...
                      'rts:Dia/RAPTOR/Phib.val','rts:Dia/RAPTOR/beta.val','rts:Dia/RAPTOR/conf_state.val', ...
                      'rts:Dia/RAPTOR/exitflag.val','rts:Dia/RAPTOR/gflag.val','rts:Dia/RAPTOR/it.val', ...
                      'rts:Dia/RAPTOR/li3.val','rts:Dia/RAPTOR/neflag.val','rts:Dia/RAPTOR/res.val', ...
                      'rts:Dia/RAPTOR/dtcomp.val','rts:Dia/RAPTOR/ECEmask.val','rts:Dia/RAPTOR/F.val', ...
                      'rts:Dia/RAPTOR/Vp.val','rts:Dia/RAPTOR/chie.val','rts:Dia/RAPTOR/dk_k.val', ...
                      'rts:Dia/RAPTOR/dpsi.val','rts:Dia/RAPTOR/dte.val','rts:Dia/RAPTOR/g1.val', ...
                      'rts:Dia/RAPTOR/g23or.val','rts:Dia/RAPTOR/g2.val','rts:Dia/RAPTOR/g3.val', ...
                      'rts:Dia/RAPTOR/iota.val','rts:Dia/RAPTOR/jbs.val','rts:Dia/RAPTOR/jec.val', ...
                      'rts:Dia/RAPTOR/jnb.val','rts:Dia/RAPTOR/jpar.val','rts:Dia/RAPTOR/ne.val', ...
                      'rts:Dia/RAPTOR/ni.val','rts:Dia/RAPTOR/p.val','rts:Dia/RAPTOR/paux.val', ...
                      'rts:Dia/RAPTOR/pbrem.val','rts:Dia/RAPTOR/pec.val','rts:Dia/RAPTOR/pei.val', ...
                      'rts:Dia/RAPTOR/pnbe.val','rts:Dia/RAPTOR/poh.val','rts:Dia/RAPTOR/prad.val', ...
                      'rts:Dia/RAPTOR/psi.val','rts:Dia/RAPTOR/q.val','rts:Dia/RAPTOR/rhotorECE.val', ...
                      'rts:Dia/RAPTOR/sdk_k.val','rts:Dia/RAPTOR/shear.val','rts:Dia/RAPTOR/signeo.val', ...
                      'rts:Dia/RAPTOR/sxk_k.val','rts:Dia/RAPTOR/te.val','rts:Dia/RAPTOR/ti.val', ...
                      'rts:Dia/RAPTOR/upl.val','rts:Dia/RAPTOR/xk_k.val','rts:Dia/RAPTOR/ze.val'};

    if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
      gdat_data.gdat_params.source = sources_avail;
    elseif ~iscell(gdat_data.gdat_params.source)
      if ischar(gdat_data.gdat_params.source)
        gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source);
        if strcmp(gdat_data.gdat_params.source(1),'o'); gdat_data.gdat_params.source = 'observer'; end
        if strcmp(gdat_data.gdat_params.source(1),'p'); gdat_data.gdat_params.source = 'predictive'; end
        if ~any(strmatch(gdat_data.gdat_params.source,lower(sources_avail)))
          if (gdat_params.nverbose>=1)
            warning(['source= ' gdat_data.gdat_params.source ' is not part of the available sources: ' sprintf('''%s'' ',sources_avail{:})]);
          end
          return
        else
          gdat_data.gdat_params.source = {gdat_data.gdat_params.source};
        end
      else
        if (gdat_params.nverbose>=1); warning([' source parameter not compatible with: ' sprintf('''%s'' ',sources_avail{:})]); end
        return
      end
    else
      for i=1:length(gdat_data.gdat_params.source)
        gdat_data.gdat_params.source{i} = lower(gdat_data.gdat_params.source{i});
        if strcmp(gdat_data.gdat_params.source{i}(1),'o'); gdat_data.gdat_params.source{i} = 'observer'; end
        if strcmp(gdat_data.gdat_params.source{i}(1),'p'); gdat_data.gdat_params.source{i} = 'predictive'; end
        if ~any(strmatch(gdat_data.gdat_params.source{i},lower(sources_avail)))
          if gdat_data.gdat_params.nverbose>=1
            warning(['source = ' gdat_data.gdat_params.source{i} ' not expected with data_request= ' data_request_eff])
          end
        end
      end
    end

    for i=1:length(gdat_data.gdat_params.source)
      it_def=0;
      for j=1:length(rap_signals)
        gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev;
        gdat_params_eff.data_request = {'RAP',[rap_signals{j} '_' gdat_data.gdat_params.source{i}(1)],gdat_params_eff.exp_name};
        abcd = gdat_aug(gdat_data.shot,gdat_params_eff);
        gdat_data.(gdat_data.gdat_params.source{i}).(rap_signals{j}) = abcd;
        gdat_data.(gdat_data.gdat_params.source{i}).(rap_signals{j}).description = rap_desc_signals{j};
        if strcmp(rap_signals{j},'q') && ~isempty(abcd.data)
          % promote q at top level of gdat_data.(gdat_data.gdat_params.source{i})
          gdat_data.(gdat_data.gdat_params.source{i}).data = abcd.data;
          gdat_data.(gdat_data.gdat_params.source{i}).t = abcd.t;
          gdat_data.(gdat_data.gdat_params.source{i}).x = abcd.x;
          gdat_data.(gdat_data.gdat_params.source{i}).units = ['q RAPTOR_' gdat_data.gdat_params.source{i}(1:3) ' various rho'];
          gdat_data.(gdat_data.gdat_params.source{i}).dim = abcd.dim;
          gdat_data.(gdat_data.gdat_params.source{i}).dimunits = abcd.dimunits;
          gdat_data.(gdat_data.gdat_params.source{i}).dimunits{1} = 'various rhotor';
          gdat_data.(gdat_data.gdat_params.source{i}).data_fullpath = abcd.data_fullpath;
          gdat_data.(gdat_data.gdat_params.source{i}).label = ['q RAPTOR_' gdat_data.gdat_params.source{i}(1:3) ' various rho'];
        end
      end
    end
    %
    % promote somethin at top level of gdat_data?
    %
    gdat_data.data_fullpath = {'RAP','signals_o in observer _p in predictive, one signal(q) copied at top',gdat_data.gdat_params.exp_name};
    % promote predictive promoted signal if non-empty, otherwise observer one
    if ~isempty(gdat_data.predictive.data)
      source_promoted = 'predictive';
    else
      source_promoted = 'observer';
    end
    gdat_data.data = gdat_data.(source_promoted).data;
    gdat_data.x = gdat_data.(source_promoted).x;
    gdat_data.t = gdat_data.(source_promoted).t;
    gdat_data.units = gdat_data.(source_promoted).units;
    gdat_data.dim = gdat_data.(source_promoted).dim;
    gdat_data.dimunits = gdat_data.(source_promoted).dimunits;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'rhotor', 'rhotor_edge', 'rhotor_norm', 'rhovol', 'volume_rho'}
    if strcmp(upper(gdat_data.gdat_params.equil),'FPG')
      gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev;
      gdat_params_eff.data_request = [{'FPG'},{'Vol'},{'AUGD'}];
      gdat_data = gdat_aug(gdat_data.shot,gdat_params_eff);
      gdat_data.gdat_params = gdat_params_prev;
      gdat_data.label = 'Vol';
      gdat_data.data_fullpath = [gdat_data.data_fullpath ' yields only edge Volume from FPG/Vol'];
      return
    end
    [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data);
    % since Lpf depends on time, need to load all first and then loop over time for easier mapping
    [phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    ndimrho = size(phi_tree.data,2);
    if ndimrho==NTIME_Lpf
      % data seems to be transposed
      ndimrho = size(phi_tree.data,1);
      itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
    else
      itotransposeback = 0;
    end
    phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
    ijnan=find(isnan(phi_tree.value));
    phi_tree.value(ijnan)=0;
    [Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
    [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
    psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
    %
    switch data_request_eff
     case {'volume_rho', 'rhovol'}
      for it=1:NTIME
        Lpf1 = Lpf1_t(it);
        psi_it=psi_tree.value(it,Lpf1:-1:1)';
        psi_axis(it)= psi_it(1);
        psi_lcfs(it)= psi_it(end);
        gdat_data.x(:,it) = sqrt(abs((psi_it-psi_axis(it)) ./(psi_lcfs(it)-psi_axis(it))));
        gdat_data.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)';
        % gdat_data.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
        gdat_data.rhovolnorm(:,it) = sqrt(abs(gdat_data.vol(:,it) ./ gdat_data.vol(end,it)));
      end
      gdat_data.dim{1} = gdat_data.x;
      gdat_data.dim{2} = gdat_data.t;
      gdat_data.dimunits{1} = 'rhopolnorm';
      gdat_data.dimunits{2} = 'time [s]';
      gdat_data.volume_edge = gdat_data.vol(end,:);
      if strcmp(data_request_eff,'volume_rho')
        gdat_data.data = gdat_data.vol;
        gdat_data.units = Vol.units;
      else strcmp(data_request_eff,'rhovol');
        gdat_data.data = gdat_data.rhovolnorm;
        gdat_data.units = ' ';
      end
     case {'rhotor', 'rhotor_edge', 'rhotor_norm'}
      b0=gdat_aug(shot,'b0');
      gdat_data.b0 = interpos(b0.t,b0.data,gdat_data.t,-1);
      for it=1:NTIME
        Lpf1 = Lpf1_t(it);
        gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)';
        gdat_data.rhotornorm(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.phi(end,it)));
        gdat_data.rhotor(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.b0(it) ./ pi));
        psi_it=psi_tree.value(it,Lpf1:-1:1)';
        psi_axis(it)= psi_it(1);
        psi_lcfs(it)= psi_it(end);
        gdat_data.x(:,it) = sqrt(abs((psi_it-psi_axis(it)) ./(psi_lcfs(it)-psi_axis(it))));
      end
      % always provide rhotor_edge so field always exists
      gdat_data.rhotor_edge = gdat_data.rhotor(end,:);
      if strcmp(data_request_eff,'rhotor')
        gdat_data.data = gdat_data.rhotor;
        gdat_data.units = 'm';
        gdat_data.dim{1} = gdat_data.x;
        gdat_data.dim{2} = gdat_data.t;
        gdat_data.dimunits{1} = 'rhopolnorm';
        gdat_data.dimunits{2} = 'time [s]';
      elseif strcmp(data_request_eff,'rhotor_edge')
        gdat_data.data = gdat_data.rhotor(end,:);
        gdat_data.units = 'm';
        gdat_data.dim{1} = gdat_data.t;
        gdat_data.dimunits{1} = 'time [s]';
      elseif strcmp(data_request_eff,'rhotor_norm')
        gdat_data.data = gdat_data.rhotornorm;
        gdat_data.units = ' ';
        gdat_data.dim{1} = gdat_data.x;
        gdat_data.dim{2} = gdat_data.t;
        gdat_data.dimunits{1} = 'rhopolnorm';
        gdat_data.dimunits{2} = 'time [s]';
      end
    end
    gdat_data.data_fullpath = [DIAG '/PFL extract in .data: ' data_request_eff];

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'sxr'}
    % sxr from B by default or else if 'source','else' is provided
    if ~isfield(gdat_data.gdat_params,'freq')|| isempty(gdat_data.gdat_params.freq)
      gdat_data.gdat_params.freq = 'slow';
    end
    if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
      gdat_data.gdat_params.source = 'G';
    elseif length(gdat_data.gdat_params.source)>1 || ...
          upper(gdat_data.gdat_params.source)<'F' || upper(gdat_data.gdat_params.source)>'M'
      if gdat_data.gdat_params.nverbose>=1
        warning(['source = ' gdat_data.gdat_params.source ' not expected with data_request= ' data_request_eff])
      end
      return
    end
    dim1_len = 'from_chord_nb'; % thus up to max(chords_ok_nb)
    if ~isfield(gdat_data.gdat_params,'camera') || isempty(gdat_data.gdat_params.camera)
      gdat_data.gdat_params.camera = [];
    elseif ischar(gdat_data.gdat_params.camera) && strcmp(gdat_data.gdat_params.camera,'central')
      gdat_data.gdat_params.source = 'J';
      gdat_data.gdat_params.camera = 49;
    elseif isnumeric(gdat_data.gdat_params.camera)
      % ok keep the array, but first dim to contain just the related chords
      dim1_len='nb_of_chords';
    else
      if gdat_data.gdat_params.nverbose>=1
        warning(['camera = ' gdat_data.gdat_params.camera ' not expected with data_request= ' data_request_eff])
      end
      return
    end
    if length(gdat_data.gdat_params.camera)==1; dim1_len='nb_of_chords'; end
    gdat_data.gdat_params.source = upper(gdat_data.gdat_params.source);
    %
    if ~isfield(gdat_data.gdat_params,'time_interval')
      gdat_data.gdat_params.time_interval = [];
    end
    [aa,bb]=unix(['ssh ' '-o "StrictHostKeyChecking no" sxaug20.aug.ipp.mpg.de WhichSX ' num2str(shot) ' ' ...
                  upper(gdat_data.gdat_params.source)]);
    chords_ok_diags=regexpi(bb,'(?<chord>[A-Z]_[0-9]+)\saddress:\s[0-9]+\sDiag:\s(?<diag>[A-Z]+)','names');
    if isempty(chords_ok_diags)
      chords_ok_nb = [];
    else
      for i=1:length(chords_ok_diags)
        chords_ok_nb(i) =  str2num(chords_ok_diags(i).chord(3:end));
      end
    end
    if isempty(gdat_data.gdat_params.camera);
      gdat_data.gdat_params.camera = chords_ok_nb;
    else
      for i=1:length(gdat_data.gdat_params.camera)
        ij=find(chords_ok_nb==gdat_data.gdat_params.camera(i));
        if ~isempty(ij)
          chords_ok_diags_tmp(i) = chords_ok_diags(ij);
        else
          % chord not in use
          chords_ok_diags_tmp(i).chord = [];
          chords_ok_diags_tmp(i).diag = [];
        end
      end
      chords_ok_diags = chords_ok_diags_tmp;
      chords_ok_nb = gdat_data.gdat_params.camera;
    end
    exp_name_eff = 'AUGD';
    icount = 0;
    nnth = 1;
    if isnumeric(gdat_data.gdat_params.freq) && gdat_data.gdat_params.freq>1;
      nnth = floor(gdat_data.gdat_params.freq+0.5);
      gdat_data.gdat_params.freq = nnth;
    end
    for i=1:length(gdat_data.gdat_params.camera)
      if ischar(gdat_data.gdat_params.freq) && strcmp(gdat_data.gdat_params.freq,'slow'); chords_ok_diags(i).diag = 'SSX'; end
      if ~isempty(chords_ok_diags(i).diag) && ~isempty(chords_ok_diags(i).chord)
        [a,e]=rdaAUG_eff(shot,chords_ok_diags(i).diag,chords_ok_diags(i).chord,exp_name_eff, ...
          gdat_data.gdat_params.time_interval,gdat_data.gdat_params.extra_arg_sf2sig);
      else
        a.data = [];
        a.t = [];
      end
      if ~isempty(a.data)
        icount = icount + 1;
        if icount == 1
          % first time has data
          % assume all chords have same time base even if from different shotfile
          % time missing one point
          if length(a.value) == length(a.t)+1
            a.t=linspace(a.range(1),a.range(2),length(a.value));
            a_time.index(2) = length(a.value);
          end
          gdat_data.t = a.t(1:nnth:end);
          gdat_data.units = a.units;
          if strcmp(dim1_len,'from_chord_nb')
            gdat_data.data = NaN*ones(max(chords_ok_nb),length(gdat_data.t));
            gdat_data.dim{1} = [1:max(chords_ok_nb)]; % simpler to have index corresponding to chord number, except for central
          else
            gdat_data.data = NaN*ones(length(chords_ok_nb),length(gdat_data.t));
            gdat_data.dim{1} = chords_ok_nb;
          end
          gdat_data.dim{2} = gdat_data.t;
          gdat_data.dimunits = [{'chord nb'}; {'s'}];
          gdat_data.data_fullpath = ['sxr from source=''' gdat_data.gdat_params.source ' uses shotfiles SX' ...
                    setdiff(unique(strvcat(chords_ok_diags(:).diag)),'SX')'];
          gdat_data.label = ['SXR/' upper(gdat_data.gdat_params.source)];
        end
        try
          if strcmp(dim1_len,'from_chord_nb')
            gdat_data.data(chords_ok_nb(i),:) = a.data(1:nnth:end);
          else
            gdat_data.data(i,:) = a.data(1:nnth:end);
          end
        catch
          if (gdat_params.nverbose>=1); disp(['problem with ' chords_ok_diags(i).diag ' ; ' chords_ok_diags(i).chord]); end

        end
      else
        % add fields not yet defined in case all cases have empty data
      end
    end
    gdat_data.chords = chords_ok_nb;
    if length(gdat_data.dim)>=1
      gdat_data.x = gdat_data.dim{1};
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'transp'}
    extra_arg_sf2sig_eff_string = '';
    if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
      extra_arg_sf2sig_eff_string = [',' gdat_data.gdat_params.extra_arg_sf2sig];
    end
    % most of the times the exp for the shotfile should be provided
    shotfile_exp_eff = gdat_params.exp_name;
    diagname='TRA';
    gdat_params.equil = 'TRE';
    gdat_data.gdat_params = gdat_params;
    TRANSP_signals;
    for i=1:size(transp_sig,1)
      if strcmp(lower(transp_sig{i,2}),'signal') || strcmp(lower(transp_sig{i,2}),'signal-group')
        try
          abcd = rdaAUG_eff(shot,diagname,transp_sig{i,1},shotfile_exp_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
          eval(['gdat_data.' transp_sig{i,1} '= abcd;']);
          % eval(['[gdat_data.' transp_sig{i,1} ',e]=rdaAUG_eff(shot,diagname,''' transp_sig{i,1} ''',shotfile_exp_eff);']);
        catch
          eval(['gdat_data.' transp_sig{i,1} '=[];']);
        end
      elseif strcmp(lower(transp_sig{i,2}),'area-base')
        clear adata_area
        try
          eval(['[adata_area]=sf2ab(diagname,shot,transp_sig{i,1},''-exp'',shotfile_exp_eff' extra_arg_sf2sig_eff_string  ');']);
          adata_area.data = adata_area.value{1};
        catch
          adata_area.value = cell(0);
        end
        eval(['gdat_data.' transp_sig{i,1} '=adata_area;']);
      elseif strcmp(lower(transp_sig{i,2}),'time-base')
        clear adata_time
        try
          eval(['[adata_time]=sf2tb(diagname,shot,transp_sig{i,1},''-exp'',shotfile_exp_eff,shotfile_exp_eff' extra_arg_sf2sig_eff_string ');']);
          adata_time.data = adata_time.value;
        catch
          adata_time.value = [];
        end
        eval(['gdat_data.' transp_sig{i,1} '=adata_time;']);
      end
    end
    % copy TIME to .t
    if isfield(gdat_data,'TIME') && isfield(gdat_data.TIME,'value') && ~isempty(gdat_data.TIME.value)
      gdat_data.t = gdat_data.TIME.value;
      gdat_data.dim{1} = gdat_data.t;
      gdat_data.dimunits{1} = gdat_data.TIME.unit;
    end

    % add equilibrium
    params_eff = gdat_params;
    params_eff.data_request = 'equil';
    gdat_data.equil = gdat_aug(shot,params_eff);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'zeff', 'z_eff'}
    zeff_max_in_main = 10.;
    zeff_min_in_main = 0.0;
    sources_avail = {'ascii','idz'};
    fields_to_copy = {'data','units','dim','dimunits','t','x','data_fullpath','label','help'};
    for i=1:length(sources_avail)
      for j=1:length(fields_to_copy)
        gdat_data.(sources_avail{i}).(fields_to_copy{j}) = [];
      end
    end
    if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
      gdat_data.gdat_params.source = sources_avail;
    elseif ~iscell(gdat_data.gdat_params.source)
      if ischar(gdat_data.gdat_params.source)
        gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source);
        if ~any(strmatch(gdat_data.gdat_params.source,sources_avail))
          if (gdat_params.nverbose>=1)
            warning(['source= ' gdat_data.gdat_params.source ' is not part of the available sources: ' sprintf('''%s'' ',sources_avail{:})]);
          end
          return
        else
          gdat_data.gdat_params.source = {gdat_data.gdat_params.source};
        end
      else
        if (gdat_params.nverbose>=1); warning([' source parameter not compatible with: ' sprintf('''%s'' ',sources_avail{:})]); end
        return
      end
    else
      for i=1:length(gdat_data.gdat_params.source)
        gdat_data.gdat_params.source{i} = lower(gdat_data.gdat_params.source{i});
        if ~any(strmatch(gdat_data.gdat_params.source{i},lower(sources_avail)))
          if gdat_data.gdat_params.nverbose>=1
            warning(['source = ' gdat_data.gdat_params.source{i} ' not expected with data_request= ' data_request_eff])
          end
        end
      end
    end
    %
    % start with ascii then IDZ, each time copy to main, so if IDZ is present it supersedes
    %
    if any(strmatch('ascii',gdat_data.gdat_params.source))
      % try ascii file from Rachael Mc Dermot in folder /afs/ipp-garching.mpg.de/home/r/rmcdermo/idl/Zeff/ascii
      % if not present ask Rachael
      dir_ascii = '/afs/ipp-garching.mpg.de/home/r/rmcdermo/idl/Zeff/ascii';
      [a,b] = unix(['ls ' dir_ascii '/*' num2str(gdat_data.shot) '*.txt']);
      if a==0
        ij1=findstr(b,'/afs');
        ij2=findstr(b,'.txt');
        file_zeff = b(ij1:ij2+3);
        aaa=textread(file_zeff,'%f');
        if mod(numel(aaa),2)~=0
          error(['problems with odd nb of points in file: ' file_zeff]);
        else
          nbpoints=length(aaa)/2;
        end
        zeff.t = aaa(1:nbpoints);
        zeff.data = aaa(nbpoints+1:end);
        gdat_data.ascii.data = min(max(zeff.data,zeff_min_in_main),zeff_max_in_main);
        gdat_data.ascii.units = '';
        gdat_data.ascii.t = zeff.t;
        gdat_data.ascii.dim{1} = gdat_data.ascii.t;
        gdat_data.ascii.dimunits{1} = 'time [s]';
        gdat_data.ascii.data_fullpath = [' from file: ' file_zeff];
        [a1,a2,a3]=fileparts(file_zeff);
        gdat_data.ascii.label = [a2 a3];
        for i=1:length(fields_to_copy)
          if isfield(gdat_data.ascii,fields_to_copy{i})
            gdat_data.(fields_to_copy{i}) = gdat_data.ascii.(fields_to_copy{i});
          end
        end
      else
        if gdat_data.gdat_params.nverbose >= 3
          disp(['No ascii file for zeff in folder: ' dir_ascii '; ask Rachael if need be'])
        end
      end
    end
    %
    % IDZ
    %
    if any(strmatch('idz',gdat_data.gdat_params.source))
      gdat_data.fit.full_path = 'PROFILES from IDZ: Zeff, unc, lower, upper';
      params_eff = gdat_params;
      params_eff.data_request={'IDZ','Zeff', params_eff.exp_name};
      try
        [gdat_data_idz_zeff,params_kin,error_status]=gdat_aug(shot,params_eff);
        if ~isempty(gdat_data_idz_zeff.data)
          gdat_data.idz.data = gdat_data_idz_zeff.data;
          gdat_data.idz.units = '';
          gdat_data.idz.t = gdat_data_idz_zeff.t;
          gdat_data.idz.x = gdat_data_idz_zeff.x;
          gdat_data.idz.dim = gdat_data_idz_zeff.dim;
          gdat_data.idz.dimunits = gdat_data_idz_zeff.dimunits;
          gdat_data.idz.dimunits{1} = 'rhopol';
          gdat_data.idz.data_fullpath = gdat_data_idz_zeff.data_fullpath;
          gdat_data.idz.label = sprintf('%s/',gdat_data_idz_zeff.data_fullpath{:});
          gdat_data.idz.label = ['mean(' gdat_data.idz.label(1:end-1) ',1)'];
          fields_to_not_copy = {'data','dim','x','dimunits'}; % only cst value at this stage to main data, should heva zeff_rho to have profile
          for i=1:length(fields_to_copy)
            if isfield(gdat_data.idz,fields_to_copy{i}) && ~any(strmatch(fields_to_copy{i},fields_to_not_copy))
              gdat_data.(fields_to_copy{i}) = gdat_data.idz.(fields_to_copy{i});
            end
          end
          gdat_data.data = min(max(mean(gdat_data_idz_zeff.data,1)',zeff_min_in_main),zeff_max_in_main);
          gdat_data.t = gdat_data_idz_zeff.t';
          gdat_data.dim{1} = gdat_data.t;
          gdat_data.dimunits{1} = gdat_data_idz_zeff.dimunits{2};
        else
          if gdat_data.gdat_params.nverbose >= 3
            disp(['No data in shotfile IDZ ; ask Rainer if need be'])
          end
        end
      catch ME
        if gdat_data.gdat_params.nverbose >= 3
          disp(['No shot file IDZ ; ask Rainer if need be'])
        end
      end
    end
    %
   otherwise
    if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_aug']); end
    error_status=901;
    return
  end

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

gdat_data.gdat_params.help = aug_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_data.mapping_for.aug = mapping_for_aug;
gdat_params = gdat_data.gdat_params;
error_status=0;

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% functions for portions called several times

function [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data);
% get basic params to be able to read results in EQ-like shotfiles
% M_Rmesh,N_Zmesh only needed for equil when 2D quantities are required
%

exp_name_eff = gdat_data.gdat_params.exp_name;
DIAG = [];
NTIME_Lpf = [];
NTIME = [];
Lpf1_t = [];
Lpf_SOL = [];
M_Rmesh = [];
N_Zmesh = [];

extra_arg_sf2sig_eff_string = '';
if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
  extra_arg_sf2sig_eff_string = [',' gdat_data.gdat_params.extra_arg_sf2sig];
end
shot=gdat_data.shot;
exp_name_eff = gdat_data.gdat_params.exp_name;
gdat_data.gdat_params.equil = upper(gdat_data.gdat_params.equil);
DIAG = gdat_data.gdat_params.equil; % 'EQI' by default at this stage, should become EQH?
% eval(['M_Rmesh_par = sf2par(DIAG,shot,''M'',''PARMV''' extra_arg_sf2sig_eff_string ');']);
[M_Rmesh_par,e]=rdaAUG_eff(shot,DIAG,'PARMV',gdat_data.gdat_params.exp_name,[], ...
          extra_arg_sf2sig_eff_string,['param:' 'M']);
M_Rmesh = M_Rmesh_par.value + 1; % nb of points
% eval(['N_Zmesh_par = sf2par(DIAG,shot,''N'',''PARMV''' extra_arg_sf2sig_eff_string ');']);
[N_Zmesh_par,e]=rdaAUG_eff(shot,DIAG,'PARMV',gdat_data.gdat_params.exp_name,[], ...
          extra_arg_sf2sig_eff_string,['param:' 'N']);
N_Zmesh = N_Zmesh_par.value + 1; % nb of points
% eval(['NTIME_par = sf2par(DIAG,shot,''NTIME'',''PARMV''' extra_arg_sf2sig_eff_string ');']);
[NTIME_par,e]=rdaAUG_eff(shot,DIAG,'PARMV',gdat_data.gdat_params.exp_name,[], ...
          extra_arg_sf2sig_eff_string,['param:' 'NTIME']);
NTIME = NTIME_par.value; % nb of points
Lpf_par = rdaAUG_eff(shot,DIAG,'Lpf',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
% since June, nb of time points in EQ results is not consistent with NTIME and time
% It seems the first NTIME points are correct, so use this explicitely
NTIME_Lpf = length(Lpf_par.value);
if ~isempty(NTIME) && ~isempty(NTIME_Lpf) && isnumeric(NTIME) && isnumeric(NTIME_Lpf) && (NTIME < NTIME_Lpf)
  if gdat_data.gdat_params.nverbose>=3; disp('WARNING: nb of times points smaller then equil results, use first NTIME points'); end
elseif ischar(NTIME) || ischar(NTIME_Lpf) || (NTIME > NTIME_Lpf)
  if gdat_data.gdat_params.nverbose >= 1
    if ischar(NTIME) || ischar(NTIME_Lpf)
      disp(['probably no data, NTIME_Lpf = ' NTIME_Lpf])
    else
      disp('ERROR: nb of times points LARGER then equil results')
      disp('this is unexpected, so stop there and ask Olivier.Sauter@epfl.ch')
    end
  end
  return
end
i_ref = 100000;
Lpf_tot = Lpf_par.value(1:NTIME); % nb of points: i_ref*nb_SOL points + nb_core
Lpf_SOL = fix(Lpf_tot/i_ref);
Lpf1_t = mod(Lpf_tot,i_ref)+1; % nb of Lpf points
if max(Lpf1_t) > 3000
  % most probably i_ref wrong
  i_ref = 10^fix(log10(max(Lpf1_t)));
  Lpf_tot = Lpf_par.value(1:NTIME); % nb of points: i_ref*nb_SOL points + nb_core
  Lpf_SOL = fix(Lpf_tot/i_ref);
  Lpf1_t = mod(Lpf_tot,i_ref)+1; % nb of Lpf points
end

[equil_time,e]=rdaAUG_eff(shot,DIAG,'time',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.t = equil_time.value(1:NTIME);
