gdat_tcv.m

    else
      gdat_data.gdat_params.source = 'liuqe';
    end
    for itime=1:length(time)
      time_eff = time(itime);
      % use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2
      [fnames_readresults]=read_results_for_chease(shot,time_eff,liuqe_version,3,[],[],[],zshift,0,1,nrz_out);
      cocos_read_results_for_chease = 17;
      if isempty(fnames_readresults)
        if gdat_params.nverbose>=1;
          warning(['could not create eqdsk file from read_results_for_chease with data_request= ' data_request_eff]);
        end
        return
      end
      ii = regexpi(fnames_readresults{4},'eqdsksigns');
      if ~isempty(ii)
        eqdskval=read_eqdsk(fnames_readresults{4},cocos_read_results_for_chease,0,[],[],1); % read_results now has relevant cocos LIUQE= 17
      else
        disp('wrong name check order of eqdsk outputs')
        return
      end
      for i=1:length(fnames_readresults)
        unix(['rm ' fnames_readresults{i}]);
      end
      %
      % run CHEASE if asked
      %
      cocos_in = 2;
      if strcmp(lower(gdat_data.gdat_params.source),'chease')
        % will give cocos_out (becoming cocos_in)=2 by default
        [fname_out,globalsvalues,namelist_struct,namelistfile_eff] = run_chease(2,eqdskval,cocos_read_results_for_chease);
        ij = strfind(fname_out,'EQDSK_COCOS_02.OUT');
        i = [];
        for i=1:length(ij)
          if ~isempty(ij(i)); break;end
        end
        eqdsk_cocos_in = read_eqdsk(fname_out{i},2,0);
      else
        %
        % transform to cocos=2 since read_results originally assumed it was cocos=2
        [eqdsk_cocos_in, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdskval,[cocos_read_results_for_chease cocos_in]);
      end
      %
      fnamefull = fullfile(['/tmp/' getenv('USER')],['EQDSK_' num2str(shot) 't' num2str(time_eff,'%.4f')]);
      % We still write COCOS=2 case, since closer to standard (in /tmp)
      if gdat_data.gdat_params.write==1
        write_eqdsk(fnamefull,eqdsk_cocos_in,cocos_in,[],[],[],1);
      end
      % Now gdat_tcv should return the convention from LIUQE which is COCOS=17, except if specified in option
      % create standard filename name from shot, time_eff (cocos will be added by write_eqdsk)
      cocos_out = 17;
      if isfield(gdat_data.gdat_params,'cocos') && ~isempty(gdat_data.gdat_params.cocos)
        cocos_out = gdat_data.gdat_params.cocos;
      else
        gdat_data.gdat_params.cocos = cocos_out;
      end
      [eqdsk_cocosout, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdsk_cocos_in,[cocos_in cocos_out]);
      % for several times, use array of structure for eqdsks,
      % cannot use it for psi(R,Z) in .data and .x since R, Z might be different at different times,
      % so project psi(R,Z) on Rmesh, Zmesh of 1st time
      if length(time) > 1
        if gdat_data.gdat_params.write==1
          gdat_data.eqdsk(itime) = write_eqdsk(fnamefull,eqdsk_cocosout,cocos_out);
        else
          gdat_data.eqdsk(itime) = eqdsk_cocosout;
        end
        if gdat_data.gdat_params.map_eqdsk_psirz==1
          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
          % do not map all psi(r,z) onto same mesh and leave .data, .dim empty (much faster)
        end
      else
        if gdat_data.gdat_params.write==1
          gdat_data.eqdsk = write_eqdsk(fnamefull,eqdsk_cocosout,cocos_out);
        else
          gdat_data.eqdsk = 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);
    if gdat_data.gdat_params.map_eqdsk_psirz==1
      gdat_data.data_fullpath=['psi(R,Z,t) on same R,Zmesh in .data and eqdsk(itime) from read_eqdsk from LIUQE' num2str(liuqe_version) ';zshift=' num2str(zshift)];
    else
      gdat_data.data_fullpath=['eqdsk(itime) from read_eqdsk from LIUQE' num2str(liuqe_version) ';zshift=' num2str(zshift)];
    end
    gdat_data.units = 'T m^2';
    gdat_data.dimunits = {'m','m','s'};
    gdat_data.request_description = ['data=psi, x=(R,Z), eqdsk contains eqdsk structure with which ' ...
                    'plot_eqdsk, write_eqdsk, read_eqdsk can be used'];

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'halpha'}
    channels = -1;
    if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels)
      channels = gdat_data.gdat_params.channels;
    end
    % '\base::pd:pd_001';

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'ids'}
    ids_empty_path = fullfile(fileparts(mfilename('fullpath')),'..','TCV_IMAS','ids_empty');

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

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

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

        ids_empty = eval(fname);
        rmpath(ids_empty_path);
      catch ME
        fprintf('Could not load empty template for %s\n',ids_top_name);
        fprintf('Available templates:\n');
        disp(dir(ids_empty_path));
        rmpath(ids_empty_path);
        rethrow(ME);
      end
    end
    if ~isfield(gdat_data.gdat_params,'error_bar') || isempty(gdat_data.gdat_params.error_bar)
      gdat_data.gdat_params.error_bar = 'delta';
    end
    try
      if ~isempty(shot)
        [ids_top,ids_top_description] = feval(['tcv_get_ids_' ids_top_name],shot,ids_empty,gdat_data.gdat_params);
        gdat_data.(ids_top_name) = ids_top;
        gdat_data.([ids_top_name '_description']) = ids_top_description;
      else
        gdat_data.(ids_top_name) = ids_empty;
        gdat_data.([ids_top_name '_description']) = ['shot empty so return default empty structure for ' ids_top_name];
      end
    catch ME_tcv_get_ids
      disp(['there is a problem with: tcv_get_ids_' ids_top_name ...
            ' , may be check if it exists in your path or test it by itself'])
      gdat_data.(ids_top_name) = ids_empty;
      gdat_data.([ids_top_name '_description']) = getReport(ME_tcv_get_ids);
      rethrow(ME_tcv_get_ids)
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'ec_data', 'aux', 'h_cd', 'nbi_data', 'ic_data', 'lh_data','ohm_data', 'bs_data'}
    % note: same time array for all at main.data level, then individual at .ec, .nbi levels
    % At this stage fill just eccd, later add nbi
    sources_avail = {'ec','nbi','ohm','bs'}; % can be set in parameter source
    % create empty structures for all, so in return one always have same substructres
    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,'trialindx') || gdat_data.gdat_params.trialindx < 0
      gdat_data.gdat_params.trialindx = [];
    end

    if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
      switch data_request_eff
       case 'ec_data'
        gdat_data.gdat_params.source = {'ec'};
       case 'ohm_data'
        gdat_data.gdat_params.source = {'ohm'};
       case 'bs_data'
        gdat_data.gdat_params.source = {'bs'};
       otherwise
        gdat_data.gdat_params.source = sources_avail;
      end
    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
    % create structure for icd sources from params and complete with defaults
    source_icd.ec = 'toray';
    source_icd.nbi = '';
    source_icd.ohm = 'ibs';
    source_icd.bs = 'ibs';
    for i=1:length(gdat_data.gdat_params.source)
      if ~isfield(gdat_data.gdat_params,['source_' gdat_data.gdat_params.source{i}])
        gdat_data.gdat_params.(['source_' gdat_data.gdat_params.source{i}]) = source_icd.(gdat_data.gdat_params.source{i});
      else
        source_icd.(gdat_data.gdat_params.source{i}) = gdat_data.gdat_params.(['source_' gdat_data.gdat_params.source{i}]);
      end
    end

    mdsopen(shot);
    field_for_main_data = 'cd_tot';
    % add each source in main.data, on ohm time array
    gdat_data.units = 'A';
    gdat_data.label=[]; % label was defined in tcv_mapping_request as char so replace to allow cells
    %
    if any(strmatch('ec',gdat_data.gdat_params.source))
      data_fullpath = '';
      ec_help = '';
      % EC
      if strcmp(lower(source_icd.ec),'toray')
        if isempty(gdat_data.gdat_params.trialindx)
          [pabs_gyro,icdtot,pow_dens,currentdrive_dens,rho_dep_pow,drho_pow,pmax,icdmax,currentdrive_dens_w2,rho_dep_icd,drho_icd] = astra_tcv_EC_exp(shot); % centralized function for toray nodes
        else
          [pabs_gyro,icdtot,pow_dens,currentdrive_dens,rho_dep_pow,drho_pow,pmax,icdmax,currentdrive_dens_w2,rho_dep_icd,drho_icd] = astra_tcv_EC_exp(shot,[],[],[],[],[],gdat_data.gdat_params.trialindx); % centralized function for toray nodes
        end
        ec_help = 'from toray icdint with extracting of effective Icd for given launcher depending on nb rays used';
        % All EC related quantities, each substructure should have at least fields data,x,t,units,dim,dimunits,label to be copied onto gdat_data
        launchers_label = {'1','2','3','4','5','6','7','8','9','tot'};
        launchers_grid = [1:10]';
        % power deposition related:
        ec_data.p_abs_plasma.data = pabs_gyro.data * 1e6;
        ec_data.p_abs_plasma.data(end+1,:) = sum(ec_data.p_abs_plasma.data,1,'omitnan');
        ec_data.p_abs_plasma.label = [strrep(pabs_gyro.comment,'MW','W') ' ; last index is total'];
        ec_data.p_abs_plasma.units = 'W';
        ec_data.p_abs_plasma.x = launchers_grid;
        ec_data.p_abs_plasma.t =pabs_gyro.tgrid;
        ec_data.p_abs_plasma.dim = {ec_data.p_abs_plasma.x, ec_data.p_abs_plasma.t};
        ec_data.p_abs_plasma.dimunits = {launchers_label, 's'};
        %
        ec_data.p_dens.data = pow_dens.data * 1e6;
        ec_data.p_dens.data(:,end+1,:) = sum(ec_data.p_dens.data,2,'omitnan');
        ec_data.p_dens.label = [strrep(pow_dens.comment,'MW','W') ' ; last index is total'];
        ec_data.p_dens.units = 'W/m^3';
        ec_data.p_dens.x = pow_dens.rgrid';
        ec_data.p_dens.rhotor_norm = ec_data.p_dens.x;
        ec_data.p_dens.t = pow_dens.tgrid;
        ec_data.p_dens.dim = {ec_data.p_dens.x, launchers_grid, ec_data.p_dens.t};
        ec_data.p_dens.dimunits = {'rhotor_norm', launchers_label, 's'};
        %
        ec_data.max_pow_dens.data = pmax.data * 1e6;
        ec_data.max_pow_dens.label = strrep(pmax.comment,'MW','W');
        ec_data.max_pow_dens.units = 'W/m^3';
        ec_data.max_pow_dens.x = [];
        ec_data.max_pow_dens.t = pmax.tgrid;
        ec_data.max_pow_dens.dim = {ec_data.max_pow_dens.t};
        ec_data.max_pow_dens.dimunits = {'s'};
        %
        ec_data.rho_max_pow_dens.data = rho_dep_pow.data * 1e6;
        ec_data.rho_max_pow_dens.label = strrep(rho_dep_pow.comment,'MW','W');
        ec_data.rho_max_pow_dens.units = 'rhotor_norm';
        ec_data.rho_max_pow_dens.x = [];
        ec_data.rho_max_pow_dens.t = rho_dep_pow.tgrid;
        ec_data.rho_max_pow_dens.dim = {ec_data.rho_max_pow_dens.t};
        ec_data.rho_max_pow_dens.dimunits = {'s'};
        %
        ec_data.width_pow_dens.data = drho_pow.data;
        ec_data.width_pow_dens.label = drho_pow.comment;
        ec_data.width_pow_dens.units = 'rhotor_norm';
        ec_data.width_pow_dens.x = [];
        ec_data.width_pow_dens.t = drho_pow.tgrid;
        ec_data.width_pow_dens.dim = {ec_data.width_pow_dens.t};
        ec_data.width_pow_dens.dimunits = {'s'};
        % current drive deposition related:
        ec_data.cd_tot.data = icdtot.data * 1e6;
        ec_data.cd_tot.data(end+1,:) = sum(ec_data.cd_tot.data,1,'omitnan');
        ec_data.cd_tot.label = [strrep(icdtot.comment,'MA','A') ' ; last index is total'];
        ec_data.cd_tot.units = 'A';
        ec_data.cd_tot.x = launchers_grid;
        ec_data.cd_tot.t = icdtot.tgrid;
        ec_data.cd_tot.dim = {ec_data.cd_tot.x, ec_data.cd_tot.t};
        ec_data.cd_tot.dimunits = {launchers_label, 's'};
        %
        ec_data.cd_dens.data = currentdrive_dens.data * 1e6;
        ec_data.cd_dens.data(:,end+1,:) = sum(ec_data.cd_dens.data,2,'omitnan');
        ec_data.cd_dens.label = [strrep(currentdrive_dens.comment,'MA','A') ' ; last index is total'];
        ec_data.cd_dens.units = 'A/m^2';
        ec_data.cd_dens.x = currentdrive_dens.rgrid';
        ec_data.cd_dens.rhotor_norm = ec_data.cd_dens.x;
        ec_data.cd_dens.t = currentdrive_dens.tgrid;
        ec_data.cd_dens.dim = {ec_data.cd_dens.x, launchers_grid, ec_data.cd_dens.t};
        ec_data.cd_dens.dimunits = {'rhotor_norm', launchers_label, 's'};
        %
        ec_data.max_cd_dens.data = icdmax.data * 1e6;
        ec_data.max_cd_dens.label = strrep(icdmax.comment,'MA','A');
        ec_data.max_cd_dens.units = 'A/m^2';
        ec_data.max_cd_dens.x = [];
        ec_data.max_cd_dens.t = icdmax.tgrid;
        ec_data.max_cd_dens.dim = {ec_data.max_cd_dens.t};
        ec_data.max_cd_dens.dimunits = {'s'};
        %
        ec_data.rho_max_cd_dens.data = rho_dep_icd.data;
        ec_data.rho_max_cd_dens.label = rho_dep_icd.comment;
        ec_data.rho_max_cd_dens.units = 'rhotor_norm';
        ec_data.rho_max_cd_dens.x = [];
        ec_data.rho_max_cd_dens.t = rho_dep_icd.tgrid;
        ec_data.rho_max_cd_dens.dim = {ec_data.rho_max_cd_dens.t};
        ec_data.rho_max_cd_dens.dimunits = {'s'};
        %
        ec_data.width_cd_dens.data = drho_icd.data;
        ec_data.width_cd_dens.label = drho_icd.comment;
        ec_data.width_cd_dens.units = 'rhotor_norm';
        ec_data.width_cd_dens.x = [];
        ec_data.width_cd_dens.t = drho_icd.tgrid;
        ec_data.width_cd_dens.dim = {ec_data.width_cd_dens.t};
        ec_data.width_cd_dens.dimunits = {'s'};
        %
        ec_data.cd_dens_doublewidth.data = currentdrive_dens_w2.data * 1e6;
        ec_data.cd_dens_doublewidth.label = [strrep(currentdrive_dens_w2.comment,'MA','A') ' ; last index is total'];
        for subfields={'x','rhotor_norm','t','dim','dimunits','units'}
          ec_data.cd_dens_doublewidth.(subfields{1}) = ec_data.cd_dens.(subfields{1});
        end
      else
        disp(['source_icd.ec = ' source_icd.ec ' not yet implemented, ask O. Sauter'])
        ec_data.p_abs_plasma = [];
        ec_data.p_abs_plasma(end+1,:) = [];
        ec_data.p_abs_plasma_label = [];
        ec_data.p_dens = [];
        ec_data.p_dens(end+1,:) = [];
        ec_data.p_dens_label = [];
        ec_data.max_pow_dens = [];
        ec_data.max_pow_dens_label = [];
        ec_data.rho_max_pow_dens = [];
        ec_data.rho_max_pow_dens_label = [];
        ec_data.width_pow_dens = [];
        ec_data.width_pow_dens_label = [];
        % current drive deposition related:
        ec_data.cd_tot = [];
        ec_data.cd_tot(end+1,:) = [];
        ec_data.cd_tot_label = [];
        ec_data.cd_dens = [];
        ec_data.cd_dens(:,end+1,:) = [];
        ec_data.cd_dens_label = [];
        ec_data.max_cd_dens = [];
        ec_data.max_cd_dens_label = [];
        ec_data.rho_max_cd_dens = [];
        ec_data.rho_max_cd_dens_label = [];
        ec_data.width_cd_dens = [];
        ec_data.width_cd_dens_label = [];
        ec_data.cd_dens_doublewidth = [];
        ec_data.cd_dens_doublewidth_label = [];
        ec_data.rho_tor_norm = [];
        ec_data.t = [];
        ec_data.launchers = [];
        gdat_data.ec.ec_data = ec_data;
        return
      end
      gdat_data.ec.ec_data = ec_data;
      if isempty(icdtot.data) || isempty(icdtot.tgrid) || ischar(icdtot.data)
        if (gdat_params.nverbose>=1)
          warning(['problems loading data for ' source_icd.ec ...
                    ' for data_request= ' data_request_eff]);
        end
      else
        % now default is icdtot, will depend on request and data_out param of some kind
        data_fullpath = ['from toray nodes using astra_tcv_EC_exp(shot), all results in .ec_data, subfield=' field_for_main_data ...
                         'in ec.data, .x, .t, .dim, .dimunits, .label, .units'];
        for subfields={'data','x','t','units','dim','dimunits','label'}
          gdat_data.ec.(subfields{1}) = gdat_data.ec.ec_data.(field_for_main_data).(subfields{1});
        end
        gdat_data.ec.data_fullpath = data_fullpath;
        gdat_data.ec.help = ec_help;
        % add to main, assume 1st one so just use this time base and x base
        % should find launcher tot index
        gdat_data.data(end+1,:) = gdat_data.ec.data(end,:);
        gdat_data.t = gdat_data.ec.t;
        if ischar(gdat_data.label); gdat_data.label = []; end;  % label was defined in tcv_mapping_request as char so replace 1st time
        gdat_data.label{end+1}=gdat_data.ec.label;
      end
    end
    %
    if any(strmatch('nb',gdat_data.gdat_params.source))
      NBH_in_TCV = 0;
      if shot >= 51641
        % NBI
        moderemote = mdsdata('data(\VSYSTEM::TCV_NBHDB_B["NBI:MODEREMOTE_FB"])');
        lcs_mode   = mdsdata('data(\VSYSTEM::TCV_NBHDB_I["NBI:LCS_MODE"])');
        % NBH in TCV equiv moderemote='ON' AND lcs_mode = 9
        NBH_in_TCV = strcmpi(strtrim(moderemote),'ON') && lcs_mode == 9;
      else
        % Nodes used in previous block only exist outside of Vista for shots after 51641
        if any(shot == [51458 51459 51460 51461 51463 51465 51470 51472 ... % 29.JAN.2016
                        51628 51629 51631 51632 51633 ...                   % 09.FEB.2016
                        51639 51640 ... 51641                               % 10.FEB.2016
                       ])
          NBH_in_TCV = 1;
        end
      end
      if NBH_in_TCV
        % should add reading from file at this stage ala summary Karpushov
      end
    end
    %
    if any(strmatch('ohm',gdat_data.gdat_params.source))
      data_fullpath = '';
      ohm_help = '';
      if strcmp(lower(source_icd.ohm),'ibs')
        ohm_data.cd_tot = gdat_tcv(shot,'\results::ibs:iohm');
        ohm_data.cd_tot.data = ohm_data.cd_tot.data';
        ohm_data.cd_tot.units = strrep(ohm_data.cd_tot.units,'kA','A');
        ohm_data.cd_dens = gdat_tcv(shot,'\results::ibs:johmav'); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R
        ohm_data.cd_dens.units = strrep(ohm_data.cd_dens.units,'kA','A');
        abc=get_grids_1d(ohm_data.cd_dens,1,1);
        ohm_data.cd_dens.rhotor_norm = abc.grids_1d.rhotornorm;
        ohm_data.cd_dens.rhopol_norm = abc.grids_1d.rhopolnorm;
        ohm_help = 'from IBS ohm related nodes, from Iohm=Ip-Icd-Ibs assuming stationary state';
      else
        disp(['source_icd.ohm = ' source_icd.ohm ' not yet implemented, ask O. Sauter'])
        for subfields={'cd_tot','cd_dens'};
          for subsubfields={'data','x','t','units','dim','dimunits','label','rhotor_norm','rhopol_norm'}
            ohm_data.(subfields{1}).(subsubfields{1}) = [];
          end
        end
      end
      gdat_data.ohm.ohm_data = ohm_data;
      data_fullpath = ['from ibs:ohmic related nodes, all results in .ohm_data, subfield=' field_for_main_data ...
                       'in ohm.data, .x, .t, .dim, .dimunits, .label, .units'];
      for subfields={'data','x','t','units','dim','dimunits','label'}
        gdat_data.ohm.(subfields{1}) = gdat_data.ohm.ohm_data.(field_for_main_data).(subfields{1});
      end
      if isempty(gdat_data.t); gdat_data.t = gdat_data.ohm.t; end
      gdat_data.ohm.data_fullpath = data_fullpath;
      gdat_data.ohm.help = ohm_help;
      gdat_data.data(end+1,:) = interpos(gdat_data.ohm.t,gdat_data.ohm.data,gdat_data.t,-0.1);
      gdat_data.label{end+1}=gdat_data.ohm.label;
    end
    %
    if any(strmatch('bs',gdat_data.gdat_params.source))
      data_fullpath = '';
      bs_help = '';
      if strcmp(lower(source_icd.bs),'ibs')
        bs_data.cd_tot = gdat_tcv(shot,'\results::ibs:ibs');
        bs_data.cd_tot.data = bs_data.cd_tot.data';
        bs_data.cd_tot.units = strrep(bs_data.cd_tot.units,'kA','A');
        bs_data.cd_dens = gdat_tcv(shot,'\results::ibs:jbsav'); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R
        bs_data.cd_dens.units = strrep(bs_data.cd_dens.units,'kA','A');
        abc=get_grids_1d(bs_data.cd_dens,1,1);
        bs_data.cd_dens.rhotor_norm = abc.grids_1d.rhotornorm;
        bs_data.cd_dens.rhopol_norm = abc.grids_1d.rhopolnorm;
        bs_help = 'from IBS bs related nodes, from Ibs=Ip-Icd-Ibs assuming stationary state';
      else
        disp(['source_icd.bs = ' source_icd.bs ' not yet implemented, ask O. Sauter'])
        for subfields={'cd_tot','cd_dens'};
          for subsubfields={'data','x','t','units','dim','dimunits','label','rhotor_norm','rhopol_norm'}
            bs_data.(subfields{1}).(subsubfields{1}) = [];
          end
        end
      end
      gdat_data.bs.bs_data = bs_data;
      data_fullpath = ['from ibs:bsic related nodes, all results in .bs_data, subfield=' field_for_main_data ...
                       'in bs.data, .x, .t, .dim, .dimunits, .label, .units'];
      for subfields={'data','x','t','units','dim','dimunits','label'}
        gdat_data.bs.(subfields{1}) = gdat_data.bs.bs_data.(field_for_main_data).(subfields{1});
      end
      if isempty(gdat_data.t); gdat_data.t = gdat_data.bs.t; end
      gdat_data.bs.data_fullpath = data_fullpath;
      gdat_data.bs.help = bs_help;
      gdat_data.data(end+1,:) = interpos(gdat_data.bs.t,gdat_data.bs.data,gdat_data.t,-0.1);
      gdat_data.label{end+1}=gdat_data.bs.label;
    end
    %
    % add all to last index of .data(:,i)
    gdat_data.data(end+1,:) = sum(gdat_data.data,1,'omitnan');
    gdat_data.x = [1:size(gdat_data.data,1)];
    gdat_data.label{end+1}='total';
    gdat_data.dim{1} = gdat_data.x;
    gdat_data.dim{2} = gdat_data.t;
    gdat_data.dimunits = {['index for each main source + total ' field_for_main_data], 's'};
    gdat_data.data_fullpath = 'see in individual source substructure';

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'mhd'}
    if isfield(gdat_data.gdat_params,'nfft') && ~isempty(gdat_data.gdat_params.nfft)
      % used in gdat_plot for spectrogram plot
    else
      gdat_data.gdat_params.nfft = 1024;
    end
    % load n=1, 2 and 3 Bdot from magnetic measurements
    n3.data = [];
    if shot< 50926
      n1=tdi('abs(mhdmode("LFS",1,1))');
      n2=tdi('abs(mhdmode("LFS",2,1))');
      n3=tdi('abs(mhdmode("LFS",3,1))');
      gdat_data.data_fullpath='abs(mhdmode("LFS",n,1));% for 1, 2, 3';
    else
      if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source)
        % gdat_data.gdat_params.source;
      else
        gdat_data.gdat_params.source = '23';
      end
      if length(gdat_data.gdat_params.source)>=2 && strcmp(gdat_data.gdat_params.source(1:2),'23')
        aaLFSz23_sect3=tdi('\atlas::DT196_MHD_001:channel_067');
        aaLFSz23_sect11=tdi('\atlas::DT196_MHD_001:channel_075');
        n1 = aaLFSz23_sect3;
        n1.data = aaLFSz23_sect3.data - aaLFSz23_sect11.data;
        n2 = aaLFSz23_sect3;
        n2.data = aaLFSz23_sect3.data + aaLFSz23_sect11.data;
        % n3=n1;
        gdat_data.data_fullpath='\atlas::DT196_MHD_001:channel_067 -+ \atlas::DT196_MHD_001:channel_075 for n=1,2, LFS_sect_3/11, z=23cm';
        if strcmp(gdat_data.gdat_params.source,'23full')
          % HFS from sec 3 and 11
          aaHFSz23_sect3=tdi('\atlas::DT196_MHD_002:channel_018');
          aaHFSz23_sect11=tdi('\atlas::DT196_MHD_002:channel_022');
          gdat_data.n1_HFS=aaHFSz23_sect3;
          gdat_data.n1_HFS.data = gdat_data.n1_HFS.data - aaHFSz23_sect11.data;
          gdat_data.n2_HFS=aaHFSz23_sect3;
          gdat_data.n2_HFS.data = gdat_data.n2_HFS.data + aaHFSz23_sect11.data;
          gdat_data.data_fullpath=['\atlas::DT196_MHD_001:channel_067 -+ \atlas::DT196_MHD_001:channel_075 for n=1,2, LFS_sect_3/11, z=23cm; _HFS' ...
                    ' same for sector HFS: MHD_002:channel_018-+MHD_002:channel_022'];
        end
      elseif strcmp(gdat_data.gdat_params.source(1),'0')
        aaLFSz0_sect3=tdi('\atlas::DT196_MHD_001:channel_083');
        aaLFSz0_sect11=tdi('\atlas::DT196_MHD_001:channel_091');
        n1 = aaLFSz0_sect3;
        n1.data = aaLFSz0_sect3.data - aaLFSz0_sect11.data;
        n2 = aaLFSz0_sect3;
        n2.data = aaLFSz0_sect3.data + aaLFSz0_sect11.data;
        % n3=n1;
        gdat_data.data_fullpath='\atlas::DT196_MHD_001:channel_083 -+ \atlas::DT196_MHD_001:channel_091 for n=1,2, LFS_sect_3/11, z=0cm';
        if strcmp(gdat_data.gdat_params.source,'0full')
          % sect 11 180deg from sec 3
          aaHFSz0_sect3=tdi('\atlas::DT196_MHD_002:channel_034');
          aaHFSz0_sect11=tdi('\atlas::DT196_MHD_002:channel_038');
          gdat_data.n1_HFS=aaHFSz0_sect11;
          gdat_data.n1_HFS.data = gdat_data.n1_HFS.data - aaHFSz0_sect11.data;
          gdat_data.n2_HFS=aaHFSz0_sect11;
          gdat_data.n2_HFS.data = gdat_data.n2_HFS.data + aaHFSz0_sect11.data;
          gdat_data.data_fullpath=['\atlas::DT196_MHD_001:channel_083 -+ \atlas::DT196_MHD_001:channel_091 for n=1,2, LFS_sect3/11, z=0cm; _HFS' ...
                    ' same for HFS: MHD_002:channel_034-+MHD_002:channel_038'];
        end
      else
        disp('should not be here in ''mhd'', ask O. Sauter')
        return
      end
    end
    if ~isempty(n1.data)
      gdat_data.data(:,1) = reshape(n1.data,length(n1.data),1);
      if length(n2.data)==length(n1.data); gdat_data.data(:,2) = reshape(n2.data,length(n2.data),1); end
      if length(n3.data)==length(n1.data); gdat_data.data(:,3) = reshape(n3.data,length(n3.data),1); end
      gdat_data.dim{1} = n1.dim{1};
      gdat_data.t = gdat_data.dim{1};
      gdat_data.dim{2} = [1; 2; 3];
      gdat_data.dimunits{1} = n1.dimunits{1};
      gdat_data.dimunits{2} = 'n number';
      if shot>= 50926
        gdat_data.dimunits{2} = 'n number, at this stage n3=n1';
        gdat_data.dimunits{2} = 'n number, at this stage n3 not computed';
      end
      gdat_data.units = 'T/s';
      gdat_data.request_description = 'delta_Bdot from magnetic probes to get n=1/odd, 2/even and 3';
      gdat_data.label = {'n=1','n=2','n=3'}; % can be used in legend(gdat_data.label)
      if shot>= 50926
        gdat_data.label = {'n odd','n even'}; % can be used in legend(gdat_data.label)
      end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   case {'ne','te'}
    % ne or Te from Thomson data on raw z mesh vs (z,t)
    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);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   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
    if liuqe_matlab==1 && strcmp(substr_liuqe,'_1'); substr_liuqe = ''; end
    % psiscatvol obtained from linear interpolation in fact so not quite ok near axis where psi is quadratic
    recompute_psiscatvol_always = 1;
    if liuqe_version==-1; recompute_psiscatvol_always = 1; end
    if all(abs(zshift)<1e-5) && liuqe_matlab==0 && recompute_psiscatvol_always==0
      psi_max=gdat_tcv([],['\results::thomson' edge_str_dot ':psi_max' substr_liuqe],'nverbose',gdat_params.nverbose);
      psiscatvol=gdat_tcv([],['\results::thomson' edge_str_dot ':psiscatvol' substr_liuqe],'nverbose',gdat_params.nverbose);
    else
      % calculate new psiscatvol
      if liuqe_matlab==0
        psitdi = tdi(['tcv_eq(''psi'',''LIUQE' substr_liuqe_tcv_eq ''')']);
        psiaxis = tdi(['tcv_eq(''psi_axis'',''LIUQE' substr_liuqe_tcv_eq ''')']);
      else
        psitdi = tdi(['tcv_eq(''psi'',''LIUQE.M' substr_liuqe_tcv_eq ''')']);
        psiaxis = tdi(['tcv_eq(''psi_axis'',''LIUQE.M' substr_liuqe_tcv_eq ''')']);
      end
      if numel(psitdi.dim)<1
        warning('problem with psitdi in gdat_tcv ')
        psitdi
        psiaxis
        return
      end
      rmesh=psitdi.dim{1};
      zmesh=psitdi.dim{2};
      zthom=mdsdata(['dim_of(\thomson' edge_str_dot ':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(gdat_data.t)
        zeffshift=interp1(gdat_data.t,zeffshift,psitdi.dim{3});
       otherwise
        if (gdat_params.nverbose>=1);
          disp(' bad time dimension for zshift')
          disp(['it should be 1 or ' num2str(length(gdat_data.t)) ' or ' num2str(length(psitdi.dim{3}))])
        end
      end
      for it=1:length(gdat_data.t)
        itpsitdi=iround_os(psitdi.dim{3},gdat_data.t(it));
        psirz=psitdi.data(:,:,itpsitdi);
        %psiscatvol0=interp2(rmesh,zmesh,psirz',0.9*ones(size(zthom)),zthom-zeffshift(itpsitdi),'spline'); % faster with interpos
        psiscatvol0=interpos2Dcartesian(rmesh,zmesh,psirz,0.9*ones(size(zthom)),zthom-zeffshift(itpsitdi),-0.1,-0.1);
        %psiscatvol0=interp2(rmesh,zmesh,psirz',0.9*ones(size(zthom)),zthom-zeffshift(itpsitdi),'linear');
        psiscatvol.data(it,:)=psiscatvol0;
        % since take closest psi(R,Z) from psitdi, should also take closest for psi_max and not interpolating
        itpsiaxis = iround_os(psiaxis.dim{1},gdat_data.t(it));
        psi_max.data(it,1) = psiaxis.data(itpsiaxis);
      end
      psiscatvol.dim{1} = gdat_data.t;
      psiscatvol.dim{2} = gdat_data.x;
    end
    if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data)
      for ir=1:length(psiscatvol.dim{2})
        rho2 = max(1.-psiscatvol.data(:,ir)./psi_max.data,0);
        rho(ir,:)= sqrt(rho2');
        % rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))';
      end
    else
      if gdat_params.nverbose>=1 && gdat_data.gdat_params.edge==0
        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;
    % add grids_1d to have rhotor, etc if cxrs_rho was asked for
    gdat_data = get_grids_1d(gdat_data,2,1,gdat_params.nverbose);

    %%%%%%%%%%%
    % 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_tcv.timedim~=2 | mapping_for_tcv.gdat_timedim~=2
        if (gdat_params.nverbose>=1);
          disp(['unexpected timedim in fit: data_request_eff= ' data_request_eff ...
                ', mapping_for_tcv.timedim= ' mapping_for_tcv.timedim ...
                ', mapping_for_tcv.gdat_timedim= ' mapping_for_tcv.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','ec','nbi','rad'}; % note should allow ech, nbh, ohmic in parameter sources
    % create empty structures for all, so in return one always have same substructres
    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
    % get ohmic power simply from vloop*Ip (minus sign for TCV), neglect dWp/dt could add it later, see chie_tcv_to_nodes
    % thus should take it from conf if present
    mdsopen(shot);
    ptot_ohm = tdi('\results::conf:ptot_ohm');
    if ~isempty(ptot_ohm.data) && ~ischar(ptot_ohm.data) && ~isempty(ptot_ohm.dim)
      gdat_data.ohm.data = ptot_ohm.data;
      gdat_data.ohm.t = ptot_ohm.dim{1};
      gdat_data.ohm.dim = ptot_ohm.dim;
      gdat_data.ohm.dimunits = ptot_ohm.dimunits;
      gdat_data.ohm.units = ptot_ohm.units;
      gdat_data.ohm.data_fullpath = '\results::conf:ptot_ohm';
      gdat_data.ohm.help = ptot_ohm.help;
    else
      params_eff = gdat_data.gdat_params;
      params_eff.data_request='ip'; % to make sure to use input params like liuqe option
      ip=gdat_tcv([],params_eff); %gdat_tcv to avoid plotting in case doplot=1 if using gdat and to save time
      params_eff.data_request='vloop';
      vloop=gdat_tcv([],params_eff);
      gdat_data.ohm.t = vloop.t;
      gdat_data.ohm.dim{1} = gdat_data.t;
      gdat_data.dimunits{1} = 's';
      gdat_data.units = 'W';
      tension = -1e5;
      vloop_smooth=interpos(vloop.t,vloop.data,gdat_data.ohm.t,tension);
      ip_t = interp1(ip.t,ip.data,gdat_data.ohm.t);
      gdat_data.ohm.data = -vloop_smooth.*ip_t; % TCV has wrong sign for Vloop
      gdat_data.ohm.raw_data = -vloop.data.*ip_t;
      gdat_data.ohm.data_fullpath = 'from vloop*Ip, smoothed vloop in data, unsmoothed in raw_data';
      gdat_data.ohm.data_fullpath=['from -vloop(tens=' num2str(tension,'%.0e') ')*ip, from gdat, in .data, unsmoothed in .raw_data'];
    end
    gdat_data.ohm.x=[];
    gdat_data.ohm.label='P_{ohm}';
    %
    % add each source in main.data, on ohm time array
    gdat_data.t = linspace(gdat_data.ohm.t(1),gdat_data.ohm.t(end),floor(1e4.*(gdat_data.ohm.t(end)-gdat_data.ohm.t(1))))';
    ij=[isfinite(gdat_data.ohm.data)];
    gdat_data.data(:,1) = interpos(-21,gdat_data.ohm.t(ij),gdat_data.ohm.data(ij),gdat_data.t);
    gdat_data.dim{1} = gdat_data.t;
    gdat_data.x(1) = 1;
    gdat_data.label={'P_{ohm}'};
    gdat_data.units = 'W';
    %
    if any(strmatch('ec',gdat_data.gdat_params.source))
      % EC
      nodenameeff='\results::toray.input:p_gyro';
      tracetdi=tdi(nodenameeff);
      if isempty(tracetdi.data) || isempty(tracetdi.dim) || ischar(tracetdi.data)
        if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
      else
        gdat_data.ec.data = tracetdi.data*1e3; % at this stage p_gyro is in kW'
        gdat_data.ec.units = 'W';
        gdat_data.ec.dim=tracetdi.dim;
        gdat_data.ec.dimunits=tracetdi.dimunits;
        gdat_data.ec.t=tracetdi.dim{1};
        gdat_data.ec.x=tracetdi.dim{2};
        gdat_data.ec.data_fullpath=[nodenameeff];
        gdat_data.ec.label='P_{EC}';
        gdat_data.ec.help = tracetdi.help;
        % add to main with linear interpolation and 0 for extrapolated values
        gdat_data.data(:,end+1) = interpos(-21,gdat_data.ec.t,gdat_data.ec.data(:,end),gdat_data.t);
        gdat_data.x(end+1) = size(gdat_data.data,2);
        gdat_data.label{end+1}=gdat_data.ec.label;
      end
    end
    %
    if any(strmatch('nb',gdat_data.gdat_params.source))
      NBH_in_TCV = 0;
      if shot >= 51641
        % NBI
        moderemote = mdsdata('data(\VSYSTEM::TCV_NBHDB_B["NBI:MODEREMOTE_FB"])');
        lcs_mode   = mdsdata('data(\VSYSTEM::TCV_NBHDB_I["NBI:LCS_MODE"])');
        % NBH in TCV equiv moderemote='ON' AND lcs_mode = 9
        NBH_in_TCV = strcmpi(strtrim(moderemote),'ON') && lcs_mode == 9;
      else
        % Nodes used in previous block only exist outside of Vista for shots after 51641
        if any(shot == [51458 51459 51460 51461 51463 51465 51470 51472 ... % 29.JAN.2016
                        51628 51629 51631 51632 51633 ...                   % 09.FEB.2016
                        51639 51640 ... 51641                               % 10.FEB.2016
                       ])
          NBH_in_TCV = 1;
        end
      end
      if NBH_in_TCV
        nodenameeff = '\results::NBH:POWR_TCV';
        nbh_data_tdi = tdi(nodenameeff);
        if ~isempty(nbh_data_tdi.data) && ~ischar(nbh_data_tdi.data) && ~isempty(nbh_data_tdi.dim)
          nbi_neutral_power_tot = nbh_data_tdi.data.*1e6; % in W
          nbi_neutral_power_tot = max(nbi_neutral_power_tot,0.);
          gdat_data.nbi.data = nbi_neutral_power_tot; % at this stage p_gyro is in kW'
          gdat_data.nbi.units = 'W';