function [gdat_data,gdat_params,error_status,varargout] = gdat_jet(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','JET' (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(51994,a2);  % gives input parameters as a structure, allows to call the same for many shots
%    a4=gdat('opt1',123,'opt2',[1 2 3],'shot',55379,'data_request','Ip','opt3','aAdB'); % all in pairs
%    a5=gdat(51994,'ip'); % standard call
%    a6=gdat(51994,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (note all lowercase in output)
%
% For JET, the specific trace can be given as:
%    aa==gdat(51994,[{'ppf'},{'efit'},{'xip'}],'machine','jet','doplot',1);
%
% Example to mask some hrts channels in the fit:
%    nete=gdat(92442,'nete','machine','jet','doplot',1,'fit_mask',[25:27]);
%
% Comments for local developer:
% This gdat is just a "header routine" calling the gdat for the specific machine gdat_`machine`.m which can be called
% directly, thus which should be able to treat the same type of input arguments
%

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

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

% construct main default parameters structure
gdat_params.data_request = '';
default_machine = 'jet';
default_user = getenv('USER');

gdat_params.machine=default_machine;
gdat_params.doplot = 0;
gdat_params.nverbose = 1;
gdat_params.jet_user = default_user;

% construct list of keywords from global set of keywords and specific JET 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 JET specific to all:
if ~isempty(data_request_names.jet)
  jet_names = fieldnames(data_request_names.jet);
  for i=1:length(jet_names)
    data_request_names.all.(jet_names{i}) = data_request_names.jet.(jet_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 = [];
gdat_data.help = [];

% Treat inputs:
ivarargin_first_char = 3;
data_request_eff = '';
if nargin>=2 && ischar(data_request); data_request = lower(data_request); end

gdat_data.gdat_request = data_request_names_all; % so if return early gives list of possible request names
gdat_data.gdat_params.help = jet_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;
    gdat_params = data_request;
  else
    % since data_request is char check from nb of args if it is data_request or a start of pairs
    if mod(nargin-1,2)==0
      ivarargin_first_char = 2;
    else
      ivarargin_first_char = 3;
      data_request_eff = data_request;
    end
  end
end

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

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

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

% if it is a request_keyword copy it:
if ischar(data_request_eff) || length(data_request_eff)==1
  ij=strmatch(data_request_eff,data_request_names_all,'exact');
else
  ij=[];
end
if ~isempty(ij);
  gdat_data.gdat_request = data_request_names_all{ij};
  if isfield(data_request_names.all.(data_request_names_all{ij}),'description') && ~isempty(data_request_names.all.(data_request_names_all{ij}).description)
    % copy description of keyword
    gdat_data.request_description = data_request_names.all.(data_request_names_all{ij}).description;
  end
else
  if ~iscell(data_request_eff) && ischar(data_request_eff)
    warning(['data_request: ' data_request_eff ' is not a pre-defined data_request nor a cell array with signal to get' char(10) ...
             'data_request probably not yet defined, ask O. Sauter if needed; check data_requests defined with' char(10) ...
             'aa=gdat(''machine'',''jet'');aa.gdat_request']);
    return
  end
end

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

if ~isfield(gdat_params,'data_request') || isempty(gdat_params.data_request)
  % warning('input for ''data_request'' missing from input arguments') % might be correct, asking for list of requests
  error_status=5;
  return
end
gdat_data.gdat_params = gdat_params;

% re-assign main variables to make sure use the one in gdat_data structure
shot = gdat_data.shot;
data_request_eff = gdat_data.gdat_params.data_request;
error_status = 6; % at least reached this level

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

% Specifications on how to get the data provided in jet_requests_mapping
mapping_for_jet = jet_requests_mapping(data_request_eff);
gdat_data.label = mapping_for_jet.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 = jet_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_data.mapping_for.jet = mapping_for_jet;
gdat_params = gdat_data.gdat_params;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1st treat the simplest method: "signal"
if strcmp(mapping_for_jet.method,'signal')
  ppftype = mapping_for_jet.expression{1};
  tracename = [mapping_for_jet.expression{2} '/' mapping_for_jet.expression{3}];
  [a,x,t,d,e]=rda_jet(shot,ppftype,tracename,gdat_params.jet_user);
  if e==0
    if gdat_params.nverbose>=3; disp(['error after rda_jet in signal with data_request_eff= ' data_request_eff]); end
    return
  end
  aatmp.data = a; aatmp.t = t; aatmp.x = x;
  gdat_data.data = aatmp.data;
  gdat_data.t = aatmp.t;
  gdat_data.x = aatmp.x;
  if isempty(gdat_data.data)
    return
  end
  if mapping_for_jet.timedim<=0
    % need to guess timedim
    if prod(size(aatmp.data)) == length(aatmp.data)
      mapping_for_jet.timedim = 1;
    elseif length(size(aatmp.data))==2
      % 2 true dimensions
      if length(aatmp.t) == size(aatmp.data,1)
        mapping_for_jet.timedim = 1;
      else
        mapping_for_jet.timedim = 2;
      end
    else
      % more than 2 dimensions, find the one with same length to define timedim
      mapping_for_jet.timedim = find(size(aatmp.data)==length(aatmp.t));
      if length(mapping_for_jet.timedim); mapping_for_jet.timedim = mapping_for_jet.timedim(1); end
    end
    mapping_for_jet.gdat_timedim = mapping_for_jet.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_jet.timedim} = gdat_data.t;
    gdat_data.dimunits{mapping_for_jet.timedim} = 'time [s]';
    gdat_data.dim{setdiff([1:2],mapping_for_jet.timedim)} = gdat_data.x;
  else
    gdat_data.dim{mapping_for_jet.timedim} = gdat_data.t;
    gdat_data.dimunits{mapping_for_jet.timedim} = 'time [s]';
    nbdims = length(size(gdat_data.data));
    for i=1:nbdims
      if i~=mapping_for_jet.timedim
        ieff=i;
        if i>mapping_for_jet.timedim; ieff=i-1; end
        gdat_data.dim{i} = gdat_data.x{ieff};
      end
    end
  end
  nbdims = length(gdat_data.dim);
  if isfield(aatmp,'units')
    gdat_data.units = aatmp.units;
  elseif isfield(mapping_for_jet,'units')
    gdat_data.units =mapping_for_jet.units;
  end
  % keyboard
  if mapping_for_jet.gdat_timedim>0 && mapping_for_jet.gdat_timedim ~= mapping_for_jet.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_jet.gdat_timedim-1);
    inew=[1:mapping_for_jet.gdat_timedim-1 mapping_for_jet.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_jet.timedim-1) 'it,' ...
                 repmat(':,',1,nbdims-mapping_for_jet.timedim-1) ':)'];
    dimstr_new=['(' repmat(':,',1,mapping_for_jet.gdat_timedim-1) 'it,' ...
                repmat(':,',1,nbdims-mapping_for_jet.gdat_timedim-1) ':)'];
    % eval gdat_data.data(;,:,...,it,...) = aatmp.data(:,:,:,it,...);
    for it=1:size(aatmp.data,mapping_for_jet.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_jet.gdat_timedim = mapping_for_jet.timedim;
  end
  gdat_data.data_fullpath=mapping_for_jet.expression;
  % end of method "signal"

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_jet.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_jet.expression ';'];
  eval([mapping_for_jet.expression ';']);
  if isempty(gdat_tmp) || (~isstruct(gdat_tmp) & ~isobject(gdat_tmp))
    if (gdat_params.nverbose>=1); warning(['expression does not create a gdat_tmp structure: ' mapping_for_jet.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_jet.timedim==-1;
      mapping_for_jet.timedim = nbdims;
      if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_jet.timedim = nbdims-1; end
    end
    dim_nontim = setdiff([1:nbdims],mapping_for_jet.timedim);
    ijt=find(strcmp(tmp_fieldnames,'t')==1);
    if isempty(ijt)
      gdat_data.t = gdat_data.dim{mapping_for_jet.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_jet.expression;
  end
  % end of method "expression"

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_jet.method,'switchcase')
  switch data_request_eff % not lower(...) since data_request_eff should be lower case already at this stage
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % First the request names valid for "all" machines:
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'cxrs'}
      %not yet finished, just started
      % load typical data from cxrs, Ti, ni, vtori and vpoli (if available), as well as zeff from cxrs
      % if 'fit' option is added: 'fit',1, then the fitted profiles are returned
      %
      % sub_nodes names from CXRS_get_profiles function, lower case when used in gdat
      sub_nodes = {'Ti','vTor','vPol','nC','Zeff'}; % first node is also copied into data, choose "default' one
      sub_nodes_out = lower({'Ti','vTor','vPol','ni','Zeff'}); %
      sub_nodes_units = {'eV','m/s','m/s','m^{-3}',''}; % first node is also copied into data, choose "default' one
      % use A. Karpushov routine to get profiles and then copy the data or the fitted profiles
      aa=CXRS_get_profiles; cxrs_params = aa.param;
      cxrs_params.k_plot=0; cxrs_params.k_debug=0;
      % add params from gdat call
      params_eff = gdat_data.gdat_params;
      if isfield(params_eff,'cxrs_plot') && params_eff.cxrs_plot>0
        cxrs_plot = params_eff.cxrs_plot;
      else
        cxrs_plot = 0;
      end
      gdat_data.gdat_params.cxrs_plot = cxrs_plot;
      if isfield(params_eff,'source') && ~isempty(params_eff.source)
        source = params_eff.source;
      else
        source = [1 2 3];
      end
      gdat_data.gdat_params.source = source;
      if isfield(params_eff,'time_interval') && ~isempty(params_eff.time_interval) && length(params_eff.time_interval)>=2
        time_interval = params_eff.time_interval;
        cxrs_plot=1;
      else
        time_interval = [];
      end
      gdat_data.gdat_params.time_interval = time_interval;
      gdat_data.gdat_params.cxrs_plot = cxrs_plot;
      fit_tension_default = -100.;
      if isfield(params_eff,'fit_tension')
        fit_tension = params_eff.fit_tension;
      else
        fit_tension = fit_tension_default;
      end
      if ~isstruct(fit_tension)
        fit_tension_eff.ti = fit_tension;
        fit_tension_eff.vi = fit_tension;
        fit_tension_eff.ni = fit_tension;
        fit_tension_eff.zeff = fit_tension;
        fit_tension = fit_tension_eff;
      else
        if ~isfield(fit_tension,'ti'); fit_tension.ti = fit_tension_default; end
        if ~isfield(fit_tension,'vi'); fit_tension.vi = fit_tension_default; end
        if ~isfield(fit_tension,'ni') && ~isfield(fit_tension,'nc'); fit_tension.ni = fit_tension_default; end
        if ~isfield(fit_tension,'zeff'); fit_tension.zeff = fit_tension_default; end
      end
      gdat_data.gdat_params.fit_tension = fit_tension;
      cxrs_params.prof.Ti.taus = fit_tension.ti;
      cxrs_params.prof.vi.taus = fit_tension.vi;
      cxrs_params.prof.nc.taus = fit_tension.ni;
      cxrs_params.prof.zeff.taus = fit_tension.zeff;
      cxrs_params.k_plot = cxrs_plot;
      cxrs_profiles = CXRS_get_profiles(shot,source,time_interval,cxrs_params);
      inb_times = length(cxrs_profiles.Times);
      gdat_data.cxrs_params = cxrs_profiles.param;
      if isempty(cxrs_profiles.Times) || ~isfield(cxrs_profiles,'proffit')
        if (gdat_params.nverbose>=1); warning(['problems loading data with CXRS_get_profiles for data_request= ' data_request_eff]); end
        for i=1:length(sub_nodes)
          sub_eff_out = sub_nodes_out{i};
          gdat_data.(sub_eff_out).fit.data = [];
          gdat_data.(sub_eff_out).fit.rho = [];
          gdat_data.(sub_eff_out).fit.error_bar = [];
          gdat_data.(sub_eff_out).raw.data = [];
          gdat_data.(sub_eff_out).raw.rho = [];
          gdat_data.(sub_eff_out).raw.error_bar = [];
          gdat_data.(sub_eff_out).raw.error_bar_rho = [];
          gdat_data.(sub_eff_out).raw.cxrs_system = [];
          gdat_data.time_interval = [];
          gdat_data.(sub_eff_out).units = sub_nodes_units{i};
          if i==1
            gdat_data.data = gdat_data.(sub_eff_out).fit.data;
            gdat_data.x = gdat_data.(sub_eff_out).fit.rho;
            gdat_data.error_bar = gdat_data.(sub_eff_out).fit.error_bar;
            gdat_data.units = gdat_data.(sub_eff_out).units;
          end
        end
        return
      end
      inb_channels =120; % need to change if gets bigger!!! but easier to prefill with NaNs and use the "use" part
      for i=1:length(sub_nodes)
        sub_eff = sub_nodes{i};
        sub_eff_out = sub_nodes_out{i};
        % fits
        if isfield(cxrs_profiles.proffit,sub_eff)
          gdat_data.(sub_eff_out).fit.data = cxrs_profiles.proffit.(sub_eff);
          gdat_data.(sub_eff_out).fit.rho = cxrs_profiles.proffit.([sub_eff '_rho']);
          gdat_data.(sub_eff_out).fit.error_bar = cxrs_profiles.proffit.(['d' sub_eff]);
        else
          gdat_data.(sub_eff_out).fit.data = [];
          gdat_data.(sub_eff_out).fit.rho = [];
          gdat_data.(sub_eff_out).fit.error_bar = [];
        end
        % raw data (use all data so keep same size)
        gdat_data.(sub_eff_out).raw.data = NaN * ones(inb_channels,inb_times);
        gdat_data.(sub_eff_out).raw.rho = NaN * ones(inb_channels,inb_times);
        gdat_data.(sub_eff_out).raw.error_bar = NaN * ones(inb_channels,inb_times);
        gdat_data.(sub_eff_out).raw.error_bar_rho = NaN * ones(inb_channels,inb_times);
        gdat_data.(sub_eff_out).raw.cxrs_system = NaN * ones(inb_channels,inb_times);
        gdat_data.time_interval = [];
        for it=1:inb_times
          if isfield(cxrs_profiles,sub_eff)
            nb_raw_points = length(cxrs_profiles.(sub_eff){it}.use.y);
            gdat_data.(sub_eff_out).raw.data(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.y;
            gdat_data.(sub_eff_out).raw.rho(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.x;
            gdat_data.(sub_eff_out).raw.error_bar(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.dy;
            gdat_data.(sub_eff_out).raw.error_bar_rho(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.dx;
            gdat_data.(sub_eff_out).raw.cxrs_system(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.sys;
            gdat_data.time_interval{it} = cxrs_profiles.(sub_eff){it}.t_lim;
          end
        end
        gdat_data.(sub_eff_out).units = sub_nodes_units{i};
        if i==1
          gdat_data.data = gdat_data.(sub_eff_out).fit.data;
          gdat_data.x = gdat_data.(sub_eff_out).fit.rho;
          gdat_data.error_bar = gdat_data.(sub_eff_out).fit.error_bar;
          gdat_data.units = gdat_data.(sub_eff_out).units;
        end
      end
      gdat_data.t = cxrs_profiles.proffit.time;
      gdat_data.dim = {gdat_data.x; gdat_data.t};
      if isempty(time_interval)
        gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[1 2 3],[],cxrs_params); % with cxrs_params'];
      else
        gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[1 2 3],[' num2str(time_interval) '],cxrs_params); % with cxrs_params'];
      end
      gdat_data.dimunits{1} = '';
      gdat_data.dimunits{2} = 's';

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ece', 'eced', 'ece_rho', 'eced_rho'}
      params_eff = gdat_data.gdat_params;
      channel_index = [1 96];
      if ~isfield(params_eff,'channel_index') || isempty(params_eff.channel_index)
        params_eff.channel_index = channel_index;
      else
        channel_index = params_eff.channel_index;
      end
      channel_index_list = channel_index;
      if numel(channel_index) == 2; channel_index_list=[channel_index(1):channel_index(2)]; end
      if isfield(params_eff,'time_interval') && ~isempty(params_eff.time_interval) && length(params_eff.time_interval)>=2
        time_interval = params_eff.time_interval;
      else
        time_interval = [];
      end
      gdat_data.gdat_params.time_interval = time_interval;
      fit_tension_default = -3.;
      if isfield(params_eff,'fit_tension')
        fit_tension = params_eff.fit_tension;
      else
        fit_tension = fit_tension_default;
      end
      i = channel_index_list(1);
      aa.data = [];
      for i=channel_index_list(2:end)
        if isempty(aa.data)
          params_eff.data_request = {'ppf','kk3',num2str(i,'te%.2d')};
          aa = gdat_jet(shot,params_eff);
          if ~isempty(aa.data)
            aa_data_start(i,:) = reshape(aa.data,1,numel(aa.data));
            aa.data = aa_data_start;
            params_eff.data_request = {'ppf','kk3',num2str(i,'rc%.2d')};
            gdat_data.rc = gdat_jet(shot,params_eff);
            rc_data(i,:) = reshape(gdat_data.rc.data,1,numel(gdat_data.rc.data));
            gdat_data.rc.data = rc_data;
          end
        else
          % assume same time for all
          % mds_dir = '/home/amerle/public/mds-ssh';
          mds_dir = '/home/amerle/public/mds-ssh3';
          if exist(mds_dir,'dir') && ~strcmp(which('mdsipmex'),fullfile(mds_dir,'mdsipmex.mexa64'))
            addpath(mds_dir);
          end
          mdsconnect(['ssh://' gdat_params.jet_user '@mdsplus.jetdata.eu']);
          if exist(mds_dir,'dir') && ~strcmp(which('mdsipmex'),fullfile(mds_dir,'mdsipmex.mexa64'))
            rmpath(mds_dir);
          end
          rda_data_request = ['ppf/kk3/' num2str(i,'te%.2d')];
          bb = mdsvalue(['_bb=jet("' rda_data_request '",' num2str(shot) ')']);
          if ~isempty(bb) && numel(bb)==size(aa.data,2)
            aa.data(i,:) = bb;
          end
          rda_data_request = ['ppf/kk3/' num2str(i,'rc%.2d')];
          bb = mdsvalue(['_bb=jet("' rda_data_request '",' num2str(shot) ')']);
          if ~isempty(bb) && numel(bb)==size(gdat_data.rc.data,2)
            gdat_data.rc.data(i,:) = bb;
          end
        end
      end
      gdat_data.data = aa.data;
      gdat_data.t = aa.t;
      gdat_data.units = 'eV';
      gdat_data.dim{1} = channel_index_list;
      gdat_data.dim{2} = gdat_data.t;
      gdat_data.x = gdat_data.dim{1};
      gdat_data.dimunits{1} = 'xx';
      gdat_data.dimunits{2} = 's';
      gdat_data.data_fullpath = aa.data_fullpath;
      gdat_data.label = aa.label;
      gdat_data.help = 'using jet() call for other data channels to avoid loading t many times';
      gdat_data.rc.units = 'm';
      gdat_data.rc.dim{1} = channel_index_list;
      gdat_data.rc.dim{2} = gdat_data.rc.t;
      gdat_data.rc.x = gdat_data.rc.dim{1};
      gdat_data.rc.dimunits{1} = 'xx';
      gdat_data.rc.dimunits{2} = 's';

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'eqdsk'}
      %
      time=1.; % default time
      if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time)
        time = gdat_data.gdat_params.time;
      else
        gdat_data.gdat_params.time = time;
        if (gdat_params.nverbose>=3); disp(['"time" is expected as an option, choose default time = ' num2str(time)]); end
      end
      gdat_data.gdat_params.time = time;
      gdat_data.t = time;
      if ~isfield(gdat_data.gdat_params,'equil') || isempty(gdat_data.gdat_params.equil)
        gdat_data.gdat_params.equil = 'EFIT';
      end
      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
      if abs(zshift-1) <= 1e-6
        zshift_eff = -99;
      else
        zshift_eff = zshift;
      end
      if ~isfield(gdat_data.gdat_params,'nr') || isempty(gdat_data.gdat_params.nr)
        gdat_data.gdat_params.nr = 129;
      end
      if ~isfield(gdat_data.gdat_params,'nz') || isempty(gdat_data.gdat_params.nz)
        gdat_data.gdat_params.nz = 129;
      end
      [efitdata,eqd]=geteqdskJET(shot,time,gdat_data.gdat_params.nr,gdat_data.gdat_params.nz,[],zshift_eff, ...
                                      gdat_data.gdat_params.equil,[],[],100,[],gdat_data.gdat_params.jet_user);
      if length(time) > 1
        gdat_data.eqdsk = eqd;
        for itime=1:length(time)
          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
            % 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
            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
        end
      else
        gdat_data.eqdsk = eqd{1};
        gdat_data.data = gdat_data.eqdsk.psi;
        gdat_data.dim{1} = gdat_data.eqdsk.rmesh;
        gdat_data.dim{2} = gdat_data.eqdsk.zmesh;
      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 geteqdskJET from efit ; 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'];
      % data loaded to create eqdsks, can be used in geteqdskJET for other times to avoid to reload all
      gdat_data.efitdata = efitdata;

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'mhd'}
      params_eff = gdat_data.gdat_params;
      if ~isfield(params_eff,'nfft') || isempty(params_eff.nfft)
        params_eff.nfft = 4096;
      end
      gdat_data.gdat_params = params_eff;
      params_eff.data_request={'JPF','DA','C1M-H305'};
      gdat_tmp = gdat_jet(shot,params_eff);
      gdat_data.sig1.data = reshape(gdat_tmp.data,length(gdat_tmp.data),1 );
      gdat_data.t = gdat_tmp.t;
      gdat_data.dim{1} = gdat_data.t;
      gdat_data.dim{2} = [1 2];
      gdat_data.dimunits = {'s', 'odd/even'};
      gdat_data.x = gdat_data.dim{2};
      gdat_data.sig1.data_request = params_eff.data_request;
      % other coil 180deg apart toroidally
      params_eff.data_request = {'JPF','DA','C1M-T009'};
      gdat_tmp=gdat_jet(shot,params_eff);
      gdat_data.sig2.data = reshape(gdat_tmp.data,length(gdat_tmp.data),1);
      gdat_data.sig2.data_request=params_eff.data_request;
      gdat_data.label={'n=odd','n=even'};
      gdat_data.full_path='n=1 in data(:,1) and .sig1; .sig2';
      gdat_data.units = 'T/s';
      gdat_data.label = {'n=odd','n=even'}; % can be used in legend(gdat_data.label)
      %
      if ~isempty(gdat_data.sig1.data) && ~isempty(gdat_data.sig2.data)
        gdat_data.data(:,1) = (gdat_data.sig1.data - gdat_data.sig2.data)./2;
        gdat_data.data(:,2) = (gdat_data.sig1.data + gdat_data.sig2.data)./2;
      end

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ne','te', 'nete'}
      nodenameeff_rho = 'rho';
      if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit) || ~isnumeric(gdat_data.gdat_params.fit)
        gdat_data.gdat_params.fit = 1; % default do fit
      end
      if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source)
        if strcmp(lower(gdat_data.gdat_params.source),'chain2')
          gdat_data.gdat_params.source = 'hrt2';
        end
        if strcmp(lower(gdat_data.gdat_params.source),'hrt2')
          nodenameeff_rho = []; % profiles already on rho
          gdat_data.gdat_params.fit = 0; % no need, chain2 is already a fit on rhopol
        end
      else
        gdat_data.gdat_params.source = 'hrts';
      end
      if ~isfield(gdat_data.gdat_params,'fit_ne_edge') || isempty(gdat_data.gdat_params.fit_ne_edge) || ~isnumeric(gdat_data.gdat_params.fit_ne_edge)
        gdat_data.gdat_params.fit_ne_edge = 1e18; % default edge ne value (-1 to let free)
      end
      if ~isfield(gdat_data.gdat_params,'fit_te_edge') || isempty(gdat_data.gdat_params.fit_te_edge) || ~isnumeric(gdat_data.gdat_params.fit_te_edge)
        gdat_data.gdat_params.fit_te_edge = 50; % default edge te value (-1 to let free)
      end
      params_eff = gdat_data.gdat_params;
      params_eff.doplot = 0;
      % ne and/or Te from Thomson data on raw z mesh vs (z,t)
      data_request_effi{1} = data_request_eff;
      if strcmp(data_request_eff,'nete')
        data_request_effi{1} = 'ne'; % start with ne
        data_request_effi{2} = 'te';
      else
        data_request_effi{2} = [];
      end
      params_eff.data_request = {'ppf',params_eff.source,data_request_effi{1}};
      aa = gdat_jet(shot,params_eff);
      if isempty(aa.data) || isempty(aa.t)
        if gdat_params.nverbose>=3;
          aa
          disp(['with data_request= ' data_request_eff])
        end
        return
      end
      gdat_data.data = aa.data;
      gdat_data.dim = aa.dim;
      gdat_data.t = aa.dim{2};
      gdat_data.x = aa.dim{1};
      gdat_data.dimunits=[{'R [m]'} ; {'time [s]'}];
      gdat_data.data_fullpath=[data_request_eff ' from ' params_eff.source];
      gdat_data.(data_request_effi{1}).data = aa.data;
      gdat_data.(data_request_effi{1}).t = aa.t;
      gdat_data.(data_request_effi{1}).r = aa.x;
      if any(strcmp(fieldnames(aa),'units')) && ~isempty(aa.units)
        gdat_data.units=aa.units;
      else
        if strcmp(data_request_effi{1},'ne')
          gdat_data.units = 'm^{-3}';
        elseif strcmp(data_request_effi{1},'te')
          gdat_data.units = 'eV';
        end
      end
      params_eff.data_request = {'ppf',params_eff.source,['d' data_request_effi{1}]};
      aaerr = gdat_jet(shot,params_eff);
      gdat_data.error_bar = aaerr.data;
      gdat_data.(data_request_effi{1}).error_bar = aaerr.data;
      if ~isempty(nodenameeff_rho)
        params_eff.data_request = {'ppf',params_eff.source,'rho'};
        aarho = gdat_jet(shot,params_eff);
        gdat_data.rhopol = abs(aarho.data);
        gdat_data.(data_request_effi{1}).rhopol = abs(aarho.data); % keep rho >0
        gdat_data.(data_request_effi{1}).rhopol_sign = aarho.data; % to be able to distinguish channels
      else
        gdat_data.rhopol = gdat_data.x;
        gdat_data.(data_request_effi{1}).rhopol = gdat_data.x;
        gdat_data.dimunits{1} = 'rhopol';
      end
      if length(gdat_data.x) == numel(gdat_data.rhopol)
        % gdat_data.x is already rhopol and 1D
        gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose);
      else
        aa=gdat_data; aa.x = aa.rhopol;
        gdat_data2 = get_grids_1d(aa,2,1,gdat_params.nverbose);
        gdat_data.grids_1d = gdat_data2.grids_1d;
        clear aa gdat_data2
      end
      % load te if 'nete' was asked
      if ~isempty(data_request_effi{2})
        params_eff.data_request = {'ppf',params_eff.source,data_request_effi{2}};
        aa = gdat_jet(shot,params_eff);
        if isempty(aa.data) || isempty(aa.t)
          if gdat_params.nverbose>=3;
            aa
            disp(['with data_request= ' data_request_effi{2}])
          end
          return
        end
        gdat_data.(data_request_effi{2}).data = aa.data;
        gdat_data.(data_request_effi{2}).t = aa.t;
        gdat_data.(data_request_effi{2}).r = aa.x;
        if any(strcmp(fieldnames(aa),'units')) && ~isempty(aa.units)
          gdat_data.(data_request_effi{2}).units=aa.units;
        else
          if strcmp(data_request_effi{2},'ne')
            gdat_data.(data_request_effi{2}).units = 'm^{-3}';
          elseif strcmp(data_request_effi{2},'te')
            gdat_data.(data_request_effi{2}).units = 'eV';
          end
        end
        params_eff.data_request = {'ppf',params_eff.source,['d' data_request_effi{2}]};
        aaerr = gdat_jet(shot,params_eff);
        gdat_data.(data_request_effi{2}).error_bar = aaerr.data;
        gdat_data.(data_request_effi{2}).rhopol = gdat_data.(data_request_effi{1}).rhopol;
        gdat_data.(data_request_effi{2}).rhopol_sign = gdat_data.(data_request_effi{1}).rhopol_sign;
        % construct pressure in .data
        gdat_data.data = 1.602e-19.* gdat_data.(data_request_effi{1}).data .* gdat_data.(data_request_effi{2}).data;
      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.d_over_drhotornorm = [];
      gdat_data.fit.ne.data = [];
      gdat_data.fit.ne.d_over_drhotornorm = [];
      gdat_data.fit.raw.te.data = [];
      gdat_data.fit.raw.te.rhopolnorm = [];
      gdat_data.fit.raw.ne.data = [];
      gdat_data.fit.raw.ne.rhopolnorm = [];
      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.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 isfield(gdat_data.gdat_params,'fit_ne_min')
        fit_ne_min = gdat_data.gdat_params.fit_ne_min;
      else
        fit_ne_min = 1e17;
      end
      if isfield(gdat_data.gdat_params,'fit_te_min')
        fit_te_min = gdat_data.gdat_params.fit_te_min;
      else
        fit_te_min = 5;
      end
      gdat_data.gdat_params.fit_te_min = fit_te_min;
      fit_mask = -1; % we mask on the indices, thus -1 means no mask
      if isfield(gdat_data.gdat_params,'fit_mask') && isnumeric(gdat_data.gdat_params.fit_mask) ...
          && ~isempty(gdat_data.gdat_params.fit_mask)
        % channel index to mask
        fit_mask = gdat_data.gdat_params.fit_mask;
      else
        fit_mask = fit_mask;
      end
      gdat_data.gdat_params.fit_mask = fit_mask;
      gdat_data.fit.mask = fit_mask;
      %
      if gdat_data.gdat_params.fit==1
        % add fits
        gdat_data.fit.t = gdat_data.t;
        rhopolfit = linspace(0,1,fit_nb_rho_points)';
        gdat_data.fit.rhotornorm = NaN*ones(length(rhopolfit),length(gdat_data.fit.t));
        gdat_data.fit.rhopolnorm = rhopolfit;
        % common part between ne and te
        indices_ok=[1:size(gdat_data.data,1)]';
        indices_ok=setdiff(indices_ok,fit_mask);
        for it=1:length(gdat_data.t)
          % make rhopol->rhotor transformation for each time since equilibrium might have changed
          irhook=find(gdat_data.grids_1d.rhopolnorm(indices_ok,it)>0 & gdat_data.grids_1d.rhopolnorm(indices_ok,it)<1); % no need for ~isnan
          [rhoeff isort]=sort(gdat_data.grids_1d.rhopolnorm(indices_ok(irhook),it));
          gdat_data.fit.rhotornorm(:,it)=interpos([0; rhoeff; 1], ...
          [0; gdat_data.grids_1d.rhotornorm(indices_ok(irhook(isort)),it); 1],rhopolfit,-0.1,[2 2],[0 1]);
        end
        for i=1:2
          if strcmp(data_request_effi{i},'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.rhopolnorm = NaN*ones(size(gdat_data.te.data));
            gdat_data.fit.te.data = gdat_data.fit.rhotornorm;
            gdat_data.fit.te.d_over_drhotornorm = gdat_data.fit.rhotornorm;
          elseif strcmp(data_request_effi{i},'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.rhopolnorm = NaN*ones(size(gdat_data.ne.data));
            gdat_data.fit.ne.data =gdat_data.fit.rhotornorm;
            gdat_data.fit.ne.d_over_drhotornorm = gdat_data.fit.rhotornorm;
          end
          for it=1:length(gdat_data.t)
            if strcmp(data_request_effi{i},'te')
              idatate = find(gdat_data.te.data(indices_ok,it)>1 & gdat_data.grids_1d.rhopolnorm(indices_ok,it)<=1.05);
              idatate = indices_ok(idatate);
              if length(idatate)>0
                gdat_data.fit.raw.te.rhopolnorm(idatate,it) = gdat_data.grids_1d.rhopolnorm(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.rhopolnorm(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) = teeff(ij)./0.1; end
                %
                teeff =  [teeff(1); teeff];
                te_err_eff =  [1e4; te_err_eff];
                if gdat_data.gdat_params.fit_te_edge < 0
                  [gdat_data.fit.te.data(:,it), gdat_data.fit.te.d_over_drhopolnorm(:,it)] = ...
                      interpos(rhoeff,teeff,rhopolfit,fit_tension.te,[1 0],[0 0],te_err_eff);
                else
                  [gdat_data.fit.te.data(:,it), gdat_data.fit.te.d_over_drhopolnorm(:,it)] = ...
                      interpos(rhoeff,teeff,rhopolfit,fit_tension.te,[1 2],[0 gdat_data.gdat_params.fit_te_edge],te_err_eff);
                end
                % impose minimum value
                ij=find(gdat_data.fit.te.data(:,it)<fit_te_min);
                if ~isempty(ij); gdat_data.fit.te.data(ij,it) = fit_te_min; end
              end
            elseif strcmp(data_request_effi{i},'ne')
              idatane = find(gdat_data.ne.data(indices_ok,it)>1 & gdat_data.grids_1d.rhopolnorm(indices_ok,it)<=1.05);
              idatane = indices_ok(idatane);
              if length(idatane)>0
                gdat_data.fit.raw.ne.rhopolnorm(idatane,it) = gdat_data.grids_1d.rhopolnorm(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.rhopolnorm(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) = neeff(ij)./0.1; end
                %
                neeff =  [neeff(1); neeff];
                ne_err_eff =  [1e21; ne_err_eff];
                if gdat_data.gdat_params.fit_ne_edge < 0
                  [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.d_over_drhopolnorm(:,it)] = interpos(rhoeff,neeff,rhopolfit,fit_tension.ne,[1 0],[0 ...
                    0],ne_err_eff);
                else
                  ij_in_1 = find(rhoeff<1);
                  [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.d_over_drhopolnorm(:,it)] = ...
                      interpos([rhoeff(ij_in_1); 1.],[neeff(ij_in_1); gdat_data.gdat_params.fit_ne_edge], ...
          rhopolfit,fit_tension.ne,[1 2],[0 gdat_data.gdat_params.fit_ne_edge],[ne_err_eff(ij_in_1);ne_err_eff(1)]);
                end
                % impose minimum value
                ij=find(gdat_data.fit.ne.data(:,it)<fit_ne_min);
                if ~isempty(ij); gdat_data.fit.ne.data(ij,it) = fit_ne_min; end
              end
            end
          end
        end
      end

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ni','ti'}
      nodenameeff_rho = 'rho';
      data_type_eff = data_request_eff;
      if strcmp(data_request_eff,'ni'); data_type_eff = 'dd'; end
      params_eff = gdat_data.gdat_params;
      if isfield(params_eff,'source') && ~isempty(params_eff.source)
        if strcmp(lower(params_eff.source),'chain2')
          params_eff.source = [data_request_eff(1) 'ion'];
        end
        if strcmp(lower(params_eff.source(2:end)),'ion')
          nodenameeff_rho = []; % profiles already on rho
        end
      else
        params_eff.source = [datar_equest_eff(1) 'ion']; % only chain2 at this stage
      end
      gdat_data.gdat_params = params_eff;
      params_eff.doplot = 0;
      % ne or Te from Thomson data on raw z mesh vs (z,t)
      params_eff.data_request = {'ppf',params_eff.source,data_type_eff};
      aa = gdat_jet(shot,params_eff);
      if isempty(aa.data) || isempty(aa.t)
        if gdat_params.nverbose>=3;
          aa
          disp(['with data_request= ' data_request_eff])
        end
        return
      end
      gdat_data.data = aa.data;
      gdat_data.dim = aa.dim;
      gdat_data.t = aa.dim{2};
      gdat_data.x = aa.dim{1};
      gdat_data.dimunits=[{'R [m]'} ; {'time [s]'}];
      gdat_data.data_fullpath=[data_request_eff ' from ' params_eff.source];
      gdat_data.(data_request_eff).data = aa.data;
      gdat_data.(data_request_eff).t = aa.t;
      gdat_data.(data_request_eff).r = aa.x;
      if any(strcmp(fieldnames(aa),'units')) && ~isempty(aa.units)
        gdat_data.units=aa.units;
      else
        if strcmp(data_request_eff,'ni')
          gdat_data.units = 'm^{-3}';
        elseif strcmp(data_request_eff,'ti')
          gdat_data.units = 'eV';
        end
      end
      params_eff.data_request = {'ppf',params_eff.source,['e' data_type_eff]}; % not given for dd but ok will be empty
      aaerr = gdat_jet(shot,params_eff);
      gdat_data.error_bar = aaerr.data;
      gdat_data.(data_request_eff).error_bar = aaerr.data;
      if ~isempty(nodenameeff_rho)
        params_eff.data_request = {'ppf',params_eff.source,'rho'};
        aarho = gdat_jet(shot,params_eff);
        gdat_data.rhopol = abs(aarho.data);
        gdat_data.(data_request_eff).rhopol = abs(aarho.data); % keep rho >0
      else
        gdat_data.rhopol = gdat_data.x;
        gdat_data.(data_request_eff).rhopol = gdat_data.x;
        gdat_data.dimunits{1} = 'rhopol';
      end

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$    case {'ne_rho', 'te_rho', 'nete_rho'}
% $$$     % if nete_rho, do first ne, then Te later (so fir stuff already done)
% $$$     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
% $$$     if ~(isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
% $$$          gdat_data.gdat_params.edge>0)
% $$$       gdat_data.gdat_params.edge = 0;
% $$$     end
% $$$     [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.edge,gdat_params.nverbose);
% $$$     % construct rho mesh
% $$$     edge_str_ = '';
% $$$     edge_str_dot = '';
% $$$     if gdat_data.gdat_params.edge
% $$$       edge_str_ = '_edge';
% $$$       edge_str_dot = '.edge';
% $$$     end
% $$$     psi_max=gdat([],['\results::thomson' edge_str_dot ':psi_max' substr_liuqe]);
% $$$     psiscatvol=gdat([],['\results::thomson' edge_str_dot ':psiscatvol' substr_liuqe]);
% $$$     if abs(zshift)>1e-5
% $$$       % calculate new psiscatvol
% $$$       psitdi=tdi(['\results::psi' substr_liuqe]);
% $$$       rmesh=psitdi.dim{1};
% $$$       zmesh=psitdi.dim{2};
% $$$       zthom=mdsdata('dim_of(\thomson:te,1)');
% $$$       zeffshift=zshift;
% $$$       % set zeffshift time array same as psitdi
% $$$       switch length(zeffshift)
% $$$         case 1
% $$$           zeffshift=zeffshift * ones(size(psitdi.dim{3}));
% $$$         case length(psitdi.dim{3})
% $$$           % ok
% $$$         case length(psiscatvol.dim{1})
% $$$           zeffshift=interp1(psiscatvol.dim{1},zeffshift,psitdi.dim{3});
% $$$         otherwise
% $$$          if (gdat_params.nverbose>=1);
% $$$            disp(' bad time dimension for zshift')
% $$$            disp(['it should be 1 or ' num2str(length(psiscatvol.dim{1})) ' or ' num2str(length(psitdi.dim{3}))])
% $$$          end
% $$$       end
% $$$       for it=1:length(psiscatvol.dim{1})
% $$$         itpsitdi=iround(psitdi.dim{3},psiscatvol.dim{1}(it));
% $$$         psirz=psitdi.data(:,:,itpsitdi);
% $$$         psiscatvol0=griddata(rmesh,zmesh,psirz',0.9*ones(size(zthom)),zthom-zeffshift(itpsitdi));
% $$$         psiscatvol.data(it,:)=psiscatvol0;
% $$$       end
% $$$     end
% $$$     if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data)
% $$$       for ir=1:length(psiscatvol.dim{2})
% $$$         rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))';
% $$$       end
% $$$     else
% $$$       if gdat_params.nverbose>=1; warning(['psiscatvol empty?, no rho calculated for data_request = ' data_request_eff]); end
% $$$       rho=[];
% $$$     end
% $$$     gdat_data.dim{1}=rho;
% $$$     gdat_data.dimunits=[{'sqrt(psi_norm)'} ; {'time [s]'}];
% $$$     gdat_data.x=rho;
% $$$     %%%%%%%%%%%
% $$$     % if nete_rho, copy data as .ne, get .te and put pe=e ne Te in data:
% $$$     if strcmp(data_request_eff(1:4),'nete')
% $$$       % note, now has ne.data_raw for density without fir_to_thomson_ratio
% $$$       gdat_data.ne.data = gdat_data.data;
% $$$       gdat_data.ne.data_raw = gdat_data.data_raw;
% $$$       gdat_data.ne.error_bar = gdat_data.error_bar;
% $$$       gdat_data.ne.error_bar_raw = gdat_data.error_bar_raw;
% $$$       gdat_data.ne.firrat=gdat_data.firrat;
% $$$       gdat_data.ne.units = 'm^{-3}';
% $$$       gdat_data = rmfield(gdat_data,{'firrat','data_raw','error_bar_raw'});
% $$$       %
% $$$       nodenameeff=['\results::thomson' edge_str_dot ':te'];
% $$$       tracetdi=tdi(nodenameeff);
% $$$       nodenameeff=['\results::thomson' edge_str_dot ':te; error_bar'];
% $$$       tracestd=tdi(['\results::thomson' edge_str_dot ':te:error_bar']);
% $$$       gdat_data.te.data=tracetdi.data';
% $$$       gdat_data.te.error_bar=tracestd.data';
% $$$       gdat_data.te.units = tracetdi.units;
% $$$       gdat_data.data_fullpath=['pe=1.6e-19*ne*Te in data, .ne, .te from \results::thomson' ...
% $$$                     edge_str_dot ':ne and te and projected on rhopol\_norm'];
% $$$       gdat_data.units='N/m^2; 1.6022e-19 ne Te';
% $$$       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);
% $$$     end
% $$$     %%%%%%%%%%% add fitted profiles if 'fit'>=1
% $$$     default_fit_type = 'conf';
% $$$     if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit)
% $$$       gdat_data.gdat_params.fit = 1;
% $$$     end
% $$$     % add empty structure for fit so results is always the same with or without call to fit=1 or 0
% $$$     gdat_data.fit.data = [];
% $$$     gdat_data.fit.x = [];
% $$$     gdat_data.fit.t = [];
% $$$     gdat_data.fit.units = [];
% $$$     gdat_data.fit.data_fullpath = [];
% $$$     if strcmp(data_request_eff(1:4),'nete')
% $$$       % note gdat_data.fit.data level is for pe
% $$$       gdat_data.fit.ne.data = [];
% $$$       gdat_data.fit.ne.x = [];
% $$$       gdat_data.fit.ne.t = [];
% $$$       gdat_data.fit.ne.units = [];
% $$$       gdat_data.fit.ne.data_fullpath = [];
% $$$       gdat_data.fit.te.data = [];
% $$$       gdat_data.fit.te.x = [];
% $$$       gdat_data.fit.te.t = [];
% $$$       gdat_data.fit.te.units = [];
% $$$       gdat_data.fit.te.data_fullpath = [];
% $$$     end
% $$$     if gdat_data.gdat_params.fit>0
% $$$       % default is from proffit:avg_time
% $$$       if ~(isfield(gdat_data.gdat_params,'fit_type') && ~isempty(gdat_data.gdat_params.fit_type)) || ~any(strcmp(gdat_data.gdat_params.fit_type,{'local', 'avg', 'conf'}))
% $$$         gdat_data.gdat_params.fit_type = default_fit_type;
% $$$       end
% $$$       switch gdat_data.gdat_params.fit_type
% $$$        case 'avg'
% $$$         def_proffit = '\results::proffit.avg_time:';
% $$$        case 'local'
% $$$         def_proffit = '\results::proffit.local_time:';
% $$$        case 'conf'
% $$$         def_proffit = '\results::conf:';
% $$$        otherwise
% $$$         if (gdat_params.nverbose>=1);
% $$$           disp('should not be in switch gdat_data.gdat_params.fit_type')
% $$$           disp('ask olivier.sauter@epfl.ch')
% $$$         end
% $$$         return
% $$$       end
% $$$       if strcmp(gdat_data.gdat_params.fit_type,'conf')
% $$$         nodenameeff = [def_proffit data_request_eff(1:2)];
% $$$       else
% $$$         if strcmp(data_request_eff(1:2),'ne')
% $$$           nodenameeff = [def_proffit 'neft_abs']; % do first ne if nete asked for
% $$$         elseif strcmp(data_request_eff(1:2),'te')
% $$$           nodenameeff = [def_proffit 'teft'];
% $$$         else
% $$$          if (gdat_params.nverbose>=1);
% $$$            disp(['should not be here: data_request_eff, data_request_eff(1:2)= ',data_request_eff, data_request_eff(1:2)]);
% $$$          end
% $$$         end
% $$$       end
% $$$       if isfield(gdat_data.gdat_params,'trialindx') && ~isempty(gdat_data.gdat_params.trialindx) && ...
% $$$             gdat_data.gdat_params.trialindx>=0
% $$$         nodenameeff=[nodenameeff ':trial'];
% $$$         trialindx = gdat_data.gdat_params.trialindx;
% $$$       else
% $$$         gdat_data.gdat_params.trialindx = [];
% $$$         trialindx = [];
% $$$       end
% $$$       tracetdi=tdi(nodenameeff);
% $$$       if isempty(trialindx)
% $$$         if ~isempty(tracetdi.data) && ~isempty(tracetdi.dim) && ~ischar(tracetdi.data)
% $$$           gdat_data.fit.data = tracetdi.data;
% $$$         else
% $$$           if gdat_params.nverbose>=1
% $$$             disp([nodenameeff ' is empty, thus no fits, check hldsi(shot) and may need to run anaprofs?'])
% $$$           end
% $$$           if strcmp(data_request_eff(1:4),'nete')
% $$$             gdat_data.fit.ne.data_fullpath = [nodenameeff ' is empty'];
% $$$             gdat_data.fit.te.data_fullpath = [nodenameeff ' is empty'];
% $$$           else
% $$$             gdat_data.fit.data_fullpath = [nodenameeff ' is empty'];
% $$$           end
% $$$           return
% $$$         end
% $$$       else
% $$$         if ~isempty(tracetdi.data) && size(tracetdi.data,3)>=trialindx+1
% $$$           gdat_data.fit.data = tracetdi.data(:,:,trialindx+1);
% $$$         else
% $$$           if gdat_params.nverbose>=1
% $$$             disp([nodenameeff ' with trialindx=' num2str(trialindx) ' is empty, thus no fits, check hldsi(shot) and may need to run anaprofs?'])
% $$$           end
% $$$           gdat_data.fit.data_fullpath = [nodenameeff ' with trialindx=' num2str(trialindx) ' is empty'];
% $$$           return
% $$$         end
% $$$       end
% $$$       gdat_data.fit.x=tracetdi.dim{1};
% $$$       gdat_data.fit.t=tracetdi.dim{2};
% $$$       if mapping_for_jet.timedim~=2 | mapping_for_jet.gdat_timedim~=2
% $$$         if (gdat_params.nverbose>=1);
% $$$           disp(['unexpected timedim in fit: data_request_eff= ' data_request_eff ...
% $$$                 ', mapping_for_jet.timedim= ' mapping_for_jet.timedim ...
% $$$                 ', mapping_for_jet.gdat_timedim= ' mapping_for_jet.gdat_timedim]);
% $$$         end
% $$$       end
% $$$       if any(strcmp(fieldnames(tracetdi),'units'))
% $$$         gdat_data.fit.units=tracetdi.units;
% $$$       end
% $$$       gdat_data.fit.data_fullpath = nodenameeff;
% $$$       % do te as well if nete asked for
% $$$       if strcmp(data_request_eff(1:4),'nete')
% $$$         gdat_data.fit.ne.data = gdat_data.fit.data;
% $$$         gdat_data.fit.ne.x = gdat_data.fit.x;
% $$$         gdat_data.fit.ne.t = gdat_data.fit.t;
% $$$         gdat_data.fit.ne.units = gdat_data.fit.units;
% $$$         gdat_data.fit.ne.data_fullpath = gdat_data.fit.data_fullpath;
% $$$         if strcmp(gdat_data.gdat_params.fit_type,'conf')
% $$$           nodenameeff = [def_proffit 'te'];
% $$$         else
% $$$           nodenameeff = [def_proffit 'teft'];
% $$$         end
% $$$         if ~isempty(trialindx); nodenameeff=[nodenameeff ':trial']; end
% $$$         tracetdi=tdi(nodenameeff);
% $$$         if isempty(trialindx)
% $$$           gdat_data.fit.te.data = tracetdi.data;
% $$$         else
% $$$           if ~isempty(tracetdi.data) && size(tracetdi.data,3)>=trialindx+1
% $$$             gdat_data.fit.te.data = tracetdi.data(:,:,trialindx+1);
% $$$           else
% $$$             return
% $$$           end
% $$$         end
% $$$         gdat_data.fit.te.x = gdat_data.fit.ne.x;
% $$$         gdat_data.fit.te.t = gdat_data.fit.ne.t;
% $$$         if any(strcmp(fieldnames(tracetdi),'units'))
% $$$           gdat_data.fit.te.units=tracetdi.units;
% $$$         end
% $$$         gdat_data.fit.te.data_fullpath = [nodenameeff];
% $$$         % construct pe=1.6022e-19*ne*te
% $$$         gdat_data.fit.data = 1.6022e-19.*gdat_data.fit.ne.data .* gdat_data.fit.te.data;
% $$$         gdat_data.fit.units = 'N/m^2; 1.6022e-19 ne Te';
% $$$         gdat_data.fit.data_fullpath = [gdat_data.fit.data_fullpath ' ; ' nodenameeff ' and pe in data'];
% $$$       end
% $$$     else
% $$$       gdat_data.gdat_params.fit_type = default_fit_type;
% $$$     end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'powers'}
    % note: same time array for all main, ec, ohm, nbi, ...
    % At this stage fill just ech, later add nbi
    sources_avail = {'ohm','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', ...
                     'rad_bulk_h','rad_bulk_u','rad_bulk_avg','rad_bulk_t'};
    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;
    params_eff.source=[]; % use default for individual calls
    sources_coeff = [];

    % ohmic, use its time-base
    params_eff.data_request='p_ohmic';
    try
      ohm=gdat_jet(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;
      sources_coeff(end+1) = 1; % to be added to sum at end
    else
      if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
      return
    end
    mapping_for_jet.timedim = 1; mapping_for_jet.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('nb',gdat_data.gdat_params.source))
      % nbi
      params_eff.data_request={'ppf','nbi','ptot'};
      try
    nbi=gdat_jet(shot,params_eff);
      catch
      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}';
        sources_coeff(end+1) = 1; % to be added to sum at end
      end
    end
    %
    if any(strmatch('ic',gdat_data.gdat_params.source))
      % ic
      params_eff.data_request={'ppf','rff','ptot'};
      try
    ic=gdat_jet(shot,params_eff);
        if isempty(ic.data) || isempty(ic.dim)
          error('Unable to get ppf/rff/ptot'); % Effectively a shortcut to the next catch block
        end
      catch
        params_eff.data_request={'ppf','icrh','ptot'};
        try
          ic=gdat_jet(shot,params_eff);
        catch
        end
      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}';
        sources_coeff(end+1) = 1; % to be added to sum at end
      end
    end
    if any(strmatch('rad',gdat_data.gdat_params.source))
      % radiated power
      params_eff.data_request='p_rad';
      try
    rad=gdat_jet(shot,params_eff);
      catch
        rad.data = [];
        rad.dim = [];
      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);
    gdat_data.x(end+1) =gdat_data.x(end)+1;
    gdat_data.label{end+1}='P_{rad}';
        sources_coeff(end+1) = 0; % to not be added to sum at end
      end
    end
    % add tot power
    aa=sum(gdat_data.data.*repmat(reshape(sources_coeff,1,numel(sources_coeff)),size(gdat_data.data,1),1),2);
    gdat_data.data(:,end+1) = aa;
    % gdat_data.data(:,end+1) = sum(gdat_data.data,2);
    gdat_data.label{end+1}='P_{tot} without Prad';
    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 power in .data(:,end)';

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'p_rad', 'prad'}
    params_eff = gdat_data.gdat_params;
    if ~isfield(params_eff,'source') || isempty(params_eff.source)
      params_eff.source = 'bolo';
    end
    gdat_data.gdat_params = params_eff;
    switch params_eff.source
     case 'prad'
      params_eff.data_request = {'ppf','prad','prad'};
      rad=gdat_jet(shot,params_eff);
      gdat_data.label={'P_{rad} from prad'};
     otherwise
      % 'bolo' by default
      gdat_data.gdat_params.source = 'bolo';
      params_eff.data_request = {'ppf','bolo','topi'};
      rad=gdat_jet(shot,params_eff);
      gdat_data.label={'P_{rad} from bolo/topi'};
      params_eff.data_request = {'ppf','bolo','tobh'};
      rad_bulk1=gdat_jet(shot,params_eff);
      gdat_data.rad_bulk_h = rad_bulk1.data;
      params_eff.data_request = {'ppf','bolo','tobu'};
      rad_bulk2=gdat_jet(shot,params_eff);
      gdat_data.rad_bulk_u = rad_bulk2.data;
      if ~isempty(rad_bulk1.t) || isempty(rad_bulk2.t)
        gdat_data.rad_bulk_t = rad_bulk1.t;
      else
        gdat_data.rad_bulk_t = rad_bulk2.t;
      end
      if isempty(gdat_data.rad_bulk_t)
        gdat_data.rad_bulk_avg = [];
      elseif numel(rad_bulk1.data) == numel(rad_bulk2.data)
        gdat_data.rad_bulk_avg = 0.5.*(rad_bulk1.data+rad_bulk2.data);
      else
        if isempty(rad_bulk1.data)
          abc=interpos(rad_bulk2.t,rad_bulk2.data,gdat_data.rad_bulk_t,-0.1);
          gdat_data.rad_bulk_avg = abc;
        elseif isempty(rad_bulk2.data)
          abc=interpos(rad_bulk1.t,rad_bulk1.data,gdat_data.rad_bulk_t,-0.1);
          gdat_data.rad_bulk_avg = abc;
        else
          abc=interpos(rad_bulk2.t,rad_bulk2.data,gdat_data.rad_bulk_t,-0.1);
          gdat_data.rad_bulk_avg = 0.5.*(rad_bulk1.data+abc);
        end
      end
    end
    gdat_data.t = rad.t;
    gdat_data.dim{1} = gdat_data.t;
    gdat_data.dimunits{1} = 's';
    gdat_data.data = rad.data;
    gdat_data.units = 'W';
    gdat_data.data_fullpath = rad.gdat_params.data_request;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$    case {'psi_edge'}
% $$$     % psi at edge, 0 by construction in Liuqe, thus not given
% $$$     nodenameeff=['\results::psi_axis' substr_liuqe];
% $$$     if liuqe_version==-1
% $$$       nodenameeff=[begstr 'q_psi' substr_liuqe];
% $$$     end
% $$$     tracetdi=tdi(nodenameeff);
% $$$     if isempty(tracetdi.data) || isempty(tracetdi.dim) || ~any(~isnan(tracetdi.data)) % || ischar(tracetdi.data) (to add?)
% $$$       if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
% $$$       if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
% $$$       return
% $$$     end
% $$$     gdat_data.data = tracetdi.data.*0;
% $$$     gdat_data.dim = tracetdi.dim;
% $$$     gdat_data.t = gdat_data.dim{1};
% $$$     gdat_data.data_fullpath=[' zero '];
% $$$     gdat_data.dimunits = tracetdi.dimunits;
% $$$     gdat_data.units = tracetdi.units;
% $$$     gdat_data.request_description = '0 since LIUQE construct psi to be zero at LCFS';
% $$$
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'q_rho'}
    % q profile on psi from liuqe
    nodenameeff=[{'ppf'},{'EFTM'},{'Q'}];
    nodenameeff=[{'ppf'},{'EFIT'},{'Q'}];
    ppftype = nodenameeff{1};
    tracename = [nodenameeff{2} '/' nodenameeff{3}];
    [a,x,t,d,e]=rda_jet(shot,ppftype,tracename,gdat_params.jet_user);
    if e==0
      if gdat_params.nverbose>=3; disp(['error after rda_jet in signal with data_request_eff= ' data_request_eff]); end
      return
    end
    aatmp.data = a; aatmp.t = t; aatmp.x = sqrt(x); % psi to rhopol
    gdat_data.data = aatmp.data;
    gdat_data.t = aatmp.t;
    gdat_data.x = aatmp.x;
    if isempty(gdat_data.data)
      return
    end
    gdat_data.dim = [{gdat_data.x} {gdat_data.t}];
    gdat_data.data_fullpath=[nodenameeff ' on rhopol'];
    gdat_data.dimunits{1} = '';
    gdat_data.dimunits{2} = 's';
    gdat_data.units = '';
    gdat_data.request_description = 'q(rhopol\_norm)';
    % add grids_1d to have rhotor, etc
    gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'rhotor', 'rhotor_edge', 'rhotor_norm', 'rhovol', 'volume_rho'}
    % Phi(LCFS) = int(Bphi dSphi), can use Eq.(11) of "Tokamak coordinate conventions: COCOS" paper:
    % O. Sauter, S.Yu. Medvedev, Comput. Phys. Commun. 184 (2013) 293–302
    % since cocos=17 for LIUQE we get:
    % q = -dPhi/dpsi => Phi = - int(q*dpsi) which should always have the sign of B0
    % get rhotornorm etc
    switch data_request_eff
     case {'rhotor', 'rhotor_edge', 'rhotor_norm'}
      nodenameeff=[{'ppf'},{'EFIT'},{'FTOR'}];
      ppftype = nodenameeff{1};
      tracename = [nodenameeff{2} '/' nodenameeff{3}];
      [a,x,t,d,e]=rda_jet(shot,ppftype,tracename,gdat_params.jet_user);
      % FTOR is on psinorm: move to rhopolnorm
      x = sqrt(x);
      if e==0
        if gdat_params.nverbose>=3; disp(['error after rda_jet in signal with nodenameeff= ' nodenameeff]); end
        if gdat_params.nverbose>=3; disp(['cannot get ' data_request_eff]); end
        return
      end
      phi = a;
      phi_edge_2d = ones(length(x),1)*a(end,:);
      gdat_data.rhotornorm = sqrt(phi./phi_edge_2d);
      gdat_data.phi_tor = phi;
      params_eff = gdat_data.gdat_params;
      params_eff.data_request='b0';
      b0=gdat_jet(shot,params_eff);
      gdat_data.b0 = interpos(b0.t,b0.data,t,-0.1);
      if strcmp(data_request_eff,'rhotor_edge')
        gdat_data.data = sqrt(phi(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)));
        gdat_data.rhotor_edge = gdat_data.data; % always provide rhotor_edge so field always exists
        gdat_data.dim = {t};
        gdat_data.t = gdat_data.dim{1};
        gdat_data.data_fullpath='sqrt(Phi_edge/pi/B0) from {''ppf''},{''EFIT''},{''FTOR''}';
        gdat_data.units = 'm';
        gdat_data.dimunits{1} = 's';
      elseif strcmp(data_request_eff,'rhotor')
        gdat_data.data = sqrt(phi./(ones(length(x),1)*reshape(gdat_data.b0,1,numel(gdat_data.b0)))./pi);
        gdat_data.rhotor_edge = gdat_data.data(end,:);
        gdat_data.dim{1} = x;
        gdat_data.dim{2} = t;
        gdat_data.x = x;
        gdat_data.t = gdat_data.dim{2};
        gdat_data.data_fullpath='sqrt(phitor/pi/B0)) from {''ppf''},{''EFIT''},{''FTOR''}';
        gdat_data.units = 'm';
        gdat_data.dimunits{1} = 'rhopol\_norm';
        gdat_data.dimunits{2} = 's';
      elseif strcmp(data_request_eff,'rhotor_norm')
        gdat_data.data = gdat_data.rhotornorm;
        gdat_data.rhotor_edge = sqrt(phi(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)));
        gdat_data.dim{1} = x;
        gdat_data.dim{2} = t;
        gdat_data.x = x;
        gdat_data.t = gdat_data.dim{2};
        gdat_data.data_fullpath='sqrt(phitor/phitor_edge), rhotor_edge=sqrt(phitor/B0/pi) from {''ppf''},{''EFIT''},{''FTOR''}';
        gdat_data.units = '';
        gdat_data.dimunits{1} = 'rhopol\_norm';
        gdat_data.dimunits{2} = 's';
      end
     case {'volume_rho', 'rhovol'}
      nodenameeff=[{'ppf'},{'EFIT'},{'VJAC'}]; % dVdpsi?
      ppftype = nodenameeff{1};
      tracename = [nodenameeff{2} '/' nodenameeff{3}];
      [a,x,t,d,e]=rda_jet(shot,ppftype,tracename,gdat_params.jet_user);
      nodenameeff2=[{'ppf'},{'EFIT'},{'VOLM'}];
      ppftype2 = nodenameeff2{1};
      tracename2 = [nodenameeff2{2} '/' nodenameeff2{3}];
      [a2,x2,t2,d2,e2]=rda_jet(shot,ppftype2,tracename2,gdat_params.jet_user);
      if e==0 || e2==0
        if gdat_params.nverbose>=3;
          disp(['error after rda_jet in signal with nodenameeff= ' nodenameeff]);
          disp(['or with nodenameeff2= ' nodenameeff2]);
          disp(['cannot get ' data_request_eff]);
          return
        end
      end
      volume = a2;
      dvdpsinorm = a;
      for it=1:length(t)
        [~,~,~,volpsi(:,it)]=interpos(x,a(:,it),-0.1);
        volpsi(:,it) = volpsi(:,it) .* volume(it)./volpsi(end,it); % ensure that total volume is the same as volm
      end
      gdat_data.x = sqrt(x);
      gdat_data.t = t;
      gdat_data.dim = [{gdat_data.x}, {gdat_data.t}];
      gdat_data.dimunits = [{'rhopol'}, {'s'}];
      gdat_data.volume_edge = volpsi(end,:);
      switch  data_request_eff
       case 'volume_rho'
        gdat_data.data = volpsi;
       case 'rhovol'
        gdat_data.data = sqrt(volpsi./(ones(length(gdat_data.x),1)*reshape(volpsi(end,:),1,length(gdat_data.t))));
      end
     otherwise
      disp(['data_request_eff = ' data_request_eff ' not defined in this section']);
      return
    end


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'rhovol','volume_rho','volume'}
    % volume_rho = vol(rho); volume = vol(LCFS) = vol(rho=1);
    % rhovol = sqrt(vol(rho)/vol(rho=1));
    warning('should change to jet relevant data')
    keyboard
    nodenameeff='\results::psitbx:vol';
    if liuqe_version==-1
      nodenameeff=[begstr 'vol' substr_liuqe];
    end
    tracetdi=tdi(nodenameeff);
    if isempty(tracetdi.data) || isempty(tracetdi.dim)
      % try to run psitbxput
      psitbxput_version = 1.3;
      psitbxput(psitbxput_version,shot);
      %ishot = mdsopen(shot);
      tracetdi=tdi(nodenameeff);
      if isempty(tracetdi.data) || isempty(tracetdi.dim)
        return
      end
    end
    gdat_data.units = tracetdi.units;
    if strcmp(data_request_eff,'volume')
      gdat_data.data = tracetdi.data(end,:);
      gdat_data.dim{1} = tracetdi.dim{2};
      gdat_data.data_fullpath=['\results::psitbx:vol(end,:)'];
      gdat_data.dimunits{1} = tracetdi.dimunits{2};
      gdat_data.request_description = 'volume(LCFS)=volume(rhopol=1)';
    else
      gdat_data.data = tracetdi.data;
      gdat_data.dim = tracetdi.dim;
      gdat_data.dimunits = tracetdi.dimunits;
      if strcmp(data_request_eff,'volume_rho')
        gdat_data.data_fullpath=['\results::psitbx:vol'];
        gdat_data.request_description = 'volume(rho)';
      elseif strcmp(data_request_eff,'rhovol')
        gdat_data.volume_edge = gdat_data.data(end,:);
        gdat_data.data = sqrt(gdat_data.data./repmat(reshape(gdat_data.volume_edge,1,size(gdat_data.data,2)),size(gdat_data.data,1),1));
        gdat_data.data_fullpath='sqrt(\results::psitbx:vol/vol_edge)';
        gdat_data.request_description = 'sqrt(volume(rho)/volume(edge))';
      else
         if (gdat_params.nverbose>=1)
           disp(['should not be here in vol cases with data_request = ' data_request_eff]);
         end
        return
      end
    end
      gdat_data.t = gdat_data.dim{mapping_for_jet.gdat_timedim};

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'sxr', 'sxr_s40', 'sxr_htv'}
    if strcmp(data_request_eff,'sxr') || strcmp(data_request_eff,'sxr_s40')
      % newer shots?
      % SXRs: JPF/DA/SXR-S40xx xx 1:17
      if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
        gdat_data.gdat_params.source = 'sxr-s40';
        gdat_data.gdat_params.source_main = {'jpf','da'};
      end
    elseif strcmp(data_request_eff,'sxr_htv')
      if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
        gdat_data.gdat_params.source = 'H';
        gdat_data.gdat_params.source_main = {'ppf','sxr'};
      elseif ~strmatch(lower(gdat_data.gdat_params.source),{'h','v','t'})
        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
    else
      error(['should not get here, data_request_eff = ' data_request_eff]);
    end
    if ~isfield(gdat_data.gdat_params,'time_interval')
      gdat_data.gdat_params.time_interval = [];
    end
    switch lower(gdat_data.gdat_params.source)
     case 'h'
      camera_all = [1:2:16];
     case 't'
      camera_all = [1:3:35];
     case 'v'
      camera_all = [1:3:35];
     otherwise
      camera_all = [1:17];
    end
    if ~isfield(gdat_data.gdat_params,'camera') || isempty(gdat_data.gdat_params.camera)
      camera = camera_all;
    elseif isnumeric(gdat_data.gdat_params.camera)
      camera = gdat_data.gdat_params.camera;
    elseif ischar(gdat_data.gdat_params.camera)
      if strcmp(lower(gdat_data.gdat_params.camera),'central')
        gdat_data.gdat_params.source = 'h';
        camera = 11;
      else
        if gdat_data.gdat_params.nverbose>=3; disp(['camera = ' gdat_data.gdat_params.camera ' not implemented']); end
        camera = camera_all;
      end
    end
    gdat_data.gdat_params.camera = camera;
    source = gdat_data.gdat_params.source;
    try
      gdat_data.chord_ok = [];
      for ichord=1:length(camera)
        aa = gdat_jet(shot,{gdat_data.gdat_params.source_main{1},gdat_data.gdat_params.source_main{2}, ...
                            [source num2str(camera(ichord),'%.2d')]});
        if isnumeric(aa.data) && ~isempty(aa.data) && ~isempty(aa.dim)
          if isempty(gdat_data.data)
            gdat_data.data = NaN*ones(numel(camera_all),numel(aa.dim{gdat_data.mapping_for.jet.timedim}));
            gdat_data.dim{1} = camera_all;
            gdat_data.dim{gdat_data.mapping_for.jet.gdat_timedim} = aa.dim{gdat_data.mapping_for.jet.timedim};
          end
          gdat_data.data(ichord,:) = aa.data';
          gdat_data.chord_ok(end+1) = camera(ichord);
        end
      end
    catch
      if gdat_data.gdat_params.nverbose>=1
        warning(['problem with sxr, source = ' source])
      end
      return
    end
    gdat_data.t = gdat_data.dim{gdat_data.mapping_for.jet.gdat_timedim};
    gdat_data.x = gdat_data.dim{1};
    gdat_data.units = 'au';
    gdat_data.dimunits = {'chord', 's'};
    gdat_data.label = ['sxr/' gdat_data.gdat_params.source num2str(camera(1)) ' to ' num2str(camera(end))];
    gdat_data.legend = num2str(gdat_data.chord_ok');
    gdat_data.data_fullpath = [gdat_data.label ', see .chord_ok for the effective channels ok'];

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'vloop'}
    if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
      gdat_data.gdat_params.source = 'ppf';
    end
    params_eff = gdat_data.gdat_params;
    switch lower(params_eff.source)
     case 'ppf'
      params_eff.data_request = [{'PPF'},{'MAGN'},{'VL'}];
      params_eff.label = 'Vloop from MAGN/VL(index)';
      channel_index = 3;
      tension = -1e2;
     case 'jpf'
      params_eff.data_request = [{'JPF'},{'DA'},{'C2-VLRRU'}];
      params_eff.label = 'Vloop Upper Restraint Ring Flux Flux Loop, non-integrated';
      channel_index = [];
      tension = -1e4;
     otherwise
      if params_eff.nverbose>=1
        warning(['source = ' params_eff.source ' not expected with data_request= ' data_request_eff])
      end
      return
    end
    if ~isfield(params_eff,'channel_index') || isempty(params_eff.channel_index)
      params_eff.channel_index = channel_index;
    else
      channel_index = params_eff.channel_index;
    end
    if ~isfield(params_eff,'tension') || isempty(params_eff.tension)
      params_eff.tension = tension;
    else
      tension = params_eff.tension;
    end

    aa=gdat_jet(shot,params_eff);
    gdat_data.gdat_params = rmfield(params_eff,'data_request');
    if ~isempty(params_eff.channel_index) && length(aa.data)~=numel(aa.data)
      gdat_data.data = aa.data(params_eff.channel_index,:);
    else
      gdat_data.data = aa.data;
    end
    gdat_data.units = 'V';
    gdat_data.dim{1} = aa.t;
    gdat_data.t = gdat_data.dim{1};
    gdat_data.dimunits{1} = 's';
    % add tension
    gdat_data.data_smooth=interpos(gdat_data.t,gdat_data.data,gdat_data.gdat_params.tension);

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

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

%if ishot==shot; mdsclose; end

gdat_data.gdat_params.help = jet_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_data.mapping_for.jet = mapping_for_jet;
gdat_params = gdat_data.gdat_params;
error_status=0;

return
