loadTCVdata.m

 function [trace,error,varargout]=loadTCVdata(shot,data_type,varargin)
%
% list of data_type currently available (when [_2,_3] is added, means can add _i to get Liuqe i):
% 'Ip'[_2,_3]   =  current
% 'zmag'[_2,_3] =  vertical position of the center of the plasma (magnetic axis)
% 'rmag'[_2,_3] =  radial position of the center of the plasma
% 'rcont'[_2,_3] =  R of plama boundary vs time
% 'zcont'[_2,_3] =  Z of plama boundary vs time
% 'vol'[_2,_3] =  volume of flux surfaces
% 'rhovol'[_2,_3] =  sqrt(V(:,t)/V(edge,t)), normalised rho variable based on volume of flux surfaces
% 'qrho'[_2,_3] =  q profile on rho mesh
% 'q95'[_2,_3] =  q95 vs time
% 'kappa', 'elon'[_2,_3] =  edge elongation vs time
% 'delta', 'triang'[_2,_3] =  edge averaged triangularity vs time
% 'deltatop', 'triangtop'[_2,_3] =  edge upper (top) triangularity vs time
% 'deltabot', 'triangbot'[_2,_3] =  edge lower (bottom) triangularity vs time
% 'neint' =  line-integrated electron density [m/m^3]
% 'ne'= ne raw profile on (z,t). ADD error bars in .std
% 'te'= Te raw profile on (z,t). ADD error bars in .std
% 'nerho'[_2,_3]= ne profile on (rho=sqrt(psi),time) mesh.Note rho is a 2D array as depends on time. ADD error bars in .std
% 'terho'[_2,_3]= Te profile on (rho=sqrt(psi),time) mesh.Note rho is a 2D array as depends on time. ADD error bars in .std
% 'nerhozshift'[_2,_3]= same as nerho but allows for zshift [m] in equil given by varargin{1}
% 'terhozshift'[_2,_3]= same as terho but allows for zshift [m] in equil given by varargin{1}
% 'profnerho' =  ne smoothed or fitted , vs (rho,t) (from Thomson auto fit)
% 'profterho' =  te smoothed or fitted , vs (rho,t) (from Thomson auto fit)
% 'neft' =  ne fitted from data on rho mesh (from proffit.local_time:neft_abs)
% 'teft' =  te fitted from data on rho mesh (from proffit.local_time:teft)
% 'neftav' =  ne fitted from averaged over time data on rho mesh (from proffit.avg_time:neft_abs)
% 'teftav' =  te fitted from averaged over time data on rho mesh (from proffit.avg_time:teft)
% 'ece'  =  electron cyclotron emission
% 'sxr'  =  soft x-ray emission
% 'sxR'  =  soft x-ray emission with varargout{1} option (requires varargin{4}!)
% 'MPX'  =  soft x-ray from wire chambers
% 'IOH'  =  current in ohmic transformer coils IOH_1
% 'vloop'  =  loop voltage
%
% 'xx_2 or xx_3' for Liuqe2 or 3: same as above for xx related to equilibrium
%
% INPUT:
% shot: shot number
% data_type: type of the required data
%
% Definition of varargin depends on data_type:
%
% data_type=sxr or ece:
%                  varargin{1}:  [i1 i2] : if not empty, assumes need many chords from i1 to i2
%                  varargin{2}:  channel status 1= unread yet, 0= read
%                                (for traces with many channel, enables to load additional channels,
%                                 like SXR, ECE, etc.)
%                  varargin{3}:  zmag for varargout{1} computation
%                       (can have more inputs for AUG, but not used here)
% data_type=nerhozshift or terhozshift:
%                  varargin{1}:  zshift [m] constant or (t) : positive, moves equil up (that is thomson effective z down)
%                                                time dependent: needs same time base as psitbx:psi
%
% OUTPUT:
%
% trace.data:   data
% trace.t:      time of reference
% trace.x:      space of reference
% trace.dim:    cell array of grids, trace.dim{1}, {2}, ...
% trace.dimunits: units of dimensions
% trace.units:  units of data
%
% Additional Output arguments depending on data_type
%
% data_type=sxR, ece:
%                varargout{1}: major radius: intersection/projection of the view lines with z=zmag
% data_type=MPX:
%                varargout{1}: te be determined
%
%
% function needed: mds functions-xtomo_geometry-get_xtomo_data (Furno's routines)
%                  VsxrTCVradius
% Example:
%         [ip,error]=loadTCVdata(shot,'Ip',1);
%         [sxr,error,R]=loadTCVdata(shot,'sxR',1);
%

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

% To allow multiple ways of writing a specific keyword, use data_type_eff within this routine
data_type_eff=data_type;
i_23=0;
liuqe_ext='';
if strcmp(data_type_eff(end-1:end),'_2') | strcmp(data_type_eff(end-1:end),'_3')
  i_23=2;
  liuqe_ext=data_type_eff(end-1:end);
end
% use keyword without eventual _2 or _3 extension to check for multiple possibilities
data_type_eff_noext=data_type_eff(1:end-i_23);
if ~isempty(strmatch(data_type_eff_noext,[{'ip'} {'i_p'} {'xip'}],'exact'))
  data_type_eff_noext='Ip';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Te'} {'t_e'} {'TE'} {'T_e'}],'exact'))
  data_type_eff_noext='te';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Ne'} {'n_e'} {'NE'} {'N_e'}],'exact'))
  data_type_eff_noext='ne';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Rcont'}],'exact'))
  data_type_eff_noext='rcont';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Zcont'}],'exact'))
  data_type_eff_noext='zcont';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Terho'}],'exact'))
  data_type_eff_noext='terho';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Terhozshift'}],'exact'))
  data_type_eff_noext='terhozshift';
end
if ~isempty(strmatch(data_type_eff_noext,[{'SXR'}],'exact'))
  data_type_eff_noext='sxr';
end
if ~isempty(strmatch(data_type_eff_noext,[{'ECE'}],'exact'))
  data_type_eff_noext='ece';
end
if ~isempty(strmatch(data_type_eff_noext,[{'VOL'} {'volume'}],'exact'))
  data_type_eff_noext='vol';
end
if ~isempty(strmatch(data_type_eff_noext,[{'RHOVOL'}],'exact'))
  data_type_eff_noext='rhovol';
end
if ~isempty(strmatch(data_type_eff_noext,[{'q_95'} {'Q95'}],'exact'))
  data_type_eff_noext='q95';
end
if ~isempty(strmatch(data_type_eff_noext,[{'elongation'} {'elon'}],'exact'))
  data_type_eff_noext='kappa';
end
if ~isempty(strmatch(data_type_eff_noext,[{'triangularity'} {'triang'}],'exact'))
  data_type_eff_noext='delta';
end
if ~isempty(strmatch(data_type_eff_noext,[{'deltaup'} {'deltau'} {'triangtop'} {'triangu'} {'triangup'}],'exact'))
  data_type_eff_noext='deltatop';
end
if ~isempty(strmatch(data_type_eff_noext,[{'deltalow'} {'deltal'} {'triangbot'} {'triangl'} {'trianglow'}],'exact'))
  data_type_eff_noext='deltabot';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Rmag'}],'exact'))
  data_type_eff_noext='rmag';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Zmag'}],'exact'))
  data_type_eff_noext='zmag';
end
if ~isempty(strmatch(data_type_eff_noext,[{'MPX'}],'exact'))
  data_type_eff_noext='MPX';
end
if ~isempty(strmatch(data_type_eff_noext,[{'ioh'} {'Ioh'} {'iot'} {'IOT'}],'exact'))
  data_type_eff_noext='IOH';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Vloop'} {'Vsurf'} {'vsurf'}],'exact'))
  data_type_eff_noext='vloop';
end

% some defaults

% nodes which have _2 and _3 equivalence, related to simpletdi case
liuqe23=[{'\results::i_p'} {'\results::z_axis'} {'\results::r_axis'} {'\results::q_psi'} ...
      {'\results::beta_tor'} {'\results::beta_pol'} {'\results::q_95'} {'\results::l_i'}  {'\results::delta_95'} ...
      {'\results::kappa_95'} {'\results::r_contour'} {'\results::z_contour'} {'\results::psi_axis'} ...
      {'\results::thomson:psiscatvol'} {'\results::thomson:psi_max'} {'\results::rms_error'} {'\results::total_energy'}];

% all keywords and corresponding case to run below
TCVkeywrdall=[{'Ip'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'rhovol'} {'qrho'} {'q95'} {'kappa'} ...
      {'delta'} {'deltatop'} {'deltabot'} {'neint'} ...
      {'ne'} {'te'} {'nerho'} {'terho'} {'nerhozshift'} {'terhozshift'} {'profnerho'} {'profterho'} ...
      {'neft'} {'teft'} {'neftav'} {'teftav'}  ...
      {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'}];
TCVsig.iip=strmatch('Ip',TCVkeywrdall,'exact');
TCVsig.izmag=strmatch('zmag',TCVkeywrdall,'exact');
TCVsig.irmag=strmatch('rmag',TCVkeywrdall,'exact');
TCVsig.ircont=strmatch('rcont',TCVkeywrdall,'exact');
TCVsig.izcont=strmatch('zcont',TCVkeywrdall,'exact');
TCVsig.ivol=strmatch('vol',TCVkeywrdall,'exact');
TCVsig.irhovol=strmatch('rhovol',TCVkeywrdall,'exact');
TCVsig.iqrho=strmatch('qrho',TCVkeywrdall,'exact');
TCVsig.iq95=strmatch('q95',TCVkeywrdall,'exact');
TCVsig.ikappa=strmatch('kappa',TCVkeywrdall,'exact');
TCVsig.idelta=strmatch('delta',TCVkeywrdall,'exact');
TCVsig.ideltatop=strmatch('deltatop',TCVkeywrdall,'exact');
TCVsig.ideltabot=strmatch('deltabot',TCVkeywrdall,'exact');
TCVsig.ineint=strmatch('neint',TCVkeywrdall,'exact');
TCVsig.ine=strmatch('ne',TCVkeywrdall,'exact');
TCVsig.ite=strmatch('te',TCVkeywrdall,'exact');
TCVsig.inerho=strmatch('nerho',TCVkeywrdall,'exact');
TCVsig.iterho=strmatch('terho',TCVkeywrdall,'exact');
TCVsig.inerhozshift=strmatch('nerhozshift',TCVkeywrdall,'exact');
TCVsig.iterhozshift=strmatch('terhozshift',TCVkeywrdall,'exact');
TCVsig.iprofnerho=strmatch('profnerho',TCVkeywrdall,'exact');
TCVsig.iprofterho=strmatch('profterho',TCVkeywrdall,'exact');
TCVsig.ineft=strmatch('neft',TCVkeywrdall,'exact');
TCVsig.iteft=strmatch('teft',TCVkeywrdall,'exact');
TCVsig.ineftav=strmatch('neftav',TCVkeywrdall,'exact');
TCVsig.iteftav=strmatch('teftav',TCVkeywrdall,'exact');
TCVsig.isxr=strmatch('sxr',TCVkeywrdall,'exact');
TCVsig.isxR=strmatch('sxR',TCVkeywrdall,'exact');
TCVsig.iece=strmatch('ece',TCVkeywrdall,'exact');
TCVsig.iMPX=strmatch('MPX',TCVkeywrdall,'exact');
TCVsig.iIOH=strmatch('IOH',TCVkeywrdall,'exact');
TCVsig.ivloop=strmatch('vloop',TCVkeywrdall,'exact');

% For each keyword, specify which case to use. As most common is 'simpletdi', fill in with this and change
% only indices needed. Usually use name of case same as keyword name
TCVkeywrdcase=cell(size(TCVkeywrdall));
TCVkeywrdcase(:)={'simpletdi'};
TCVkeywrdcase(TCVsig.iqrho)=TCVkeywrdall(TCVsig.iqrho); % special as liuqe q_psi on psi
TCVkeywrdcase(TCVsig.ivol)=TCVkeywrdall(TCVsig.ivol); % special as nodes _2 or _3 not existing with psitbx
TCVkeywrdcase(TCVsig.irhovol)=TCVkeywrdall(TCVsig.irhovol); % idem vol
TCVkeywrdcase(TCVsig.ine)=TCVkeywrdall(TCVsig.ine); % special as dimensions from other nodes
TCVkeywrdcase(TCVsig.ite)=TCVkeywrdall(TCVsig.ite); % idem
TCVkeywrdcase(TCVsig.inerho)=TCVkeywrdall(TCVsig.inerho); % idem
TCVkeywrdcase(TCVsig.iterho)=TCVkeywrdall(TCVsig.iterho); % idem
TCVkeywrdcase(TCVsig.inerhozshift)=TCVkeywrdall(TCVsig.inerhozshift); % idem
TCVkeywrdcase(TCVsig.iterhozshift)=TCVkeywrdall(TCVsig.iterhozshift); % idem
TCVkeywrdcase(TCVsig.iprofnerho)=TCVkeywrdall(TCVsig.iprofnerho);
TCVkeywrdcase(TCVsig.iprofterho)=TCVkeywrdall(TCVsig.iprofterho);
TCVkeywrdcase(TCVsig.isxr)=TCVkeywrdall(TCVsig.isxr);
TCVkeywrdcase(TCVsig.isxR)=TCVkeywrdall(TCVsig.isxR);
TCVkeywrdcase(TCVsig.iece)=TCVkeywrdall(TCVsig.iece);
TCVkeywrdcase(TCVsig.iMPX)=TCVkeywrdall(TCVsig.iMPX);
TCVkeywrdcase(TCVsig.iIOH)=TCVkeywrdall(TCVsig.iIOH);
TCVkeywrdcase(TCVsig.ivloop)=TCVkeywrdall(TCVsig.ivloop);

% Information about which dimension has time, always return 2D data as (x,t) array
% as most are 1D arrays with time as first index, fill in with ones and change only those needed
TCVsigtimeindx=ones(size(TCVkeywrdall));

% For the 'simpletdi' cases, we need the full node in case gdat was called with full location directly
% for the other cases, leave this location empty
TCVsiglocation=cell(size(TCVkeywrdall));
TCVsiglocation(:)={''};
TCVsiglocation(TCVsig.iip)={'\results::i_p'};
TCVsiglocation(TCVsig.izmag)={'\results::z_axis'};
TCVsiglocation(TCVsig.irmag)={'\results::r_axis'};
TCVsiglocation(TCVsig.ircont)={'\results::r_contour'}; TCVsigtimeindx(TCVsig.ircont)=2;
TCVsiglocation(TCVsig.izcont)={'\results::z_contour'}; TCVsigtimeindx(TCVsig.izcont)=2;
TCVsiglocation(TCVsig.iq95)={'\results::q_95'};
TCVsiglocation(TCVsig.ikappa)={'\results::kappa_edge'};
TCVsiglocation(TCVsig.idelta)={'\results::delta_edge'};
TCVsiglocation(TCVsig.ideltatop)={'\results::delta_ed_top'};
TCVsiglocation(TCVsig.ideltabot)={'\results::delta_ed_bot'};
TCVsiglocation(TCVsig.ineint)={'\results::fir:lin_int_dens'};
%TCVsiglocation(TCVsig.iprofnerho)={'\results::THOMSON.PROFILES.AUTO:ne'};
%TCVsiglocation(TCVsig.iprofterho)={'\results::THOMSON.PROFILES.AUTO:te'};
TCVsiglocation(TCVsig.ineft)={'\results::proffit.local_time:neft_abs'}; TCVsigtimeindx(TCVsig.ineft)=2;
TCVsiglocation(TCVsig.iteft)={'\results::proffit.local_time:teft'}; TCVsigtimeindx(TCVsig.iteft)=2;
TCVsiglocation(TCVsig.ineftav)={'\results::proffit.avg_time:neft_abs'}; TCVsigtimeindx(TCVsig.ineftav)=2;
TCVsiglocation(TCVsig.iteftav)={'\results::proffit.avg_time:teft'}; TCVsigtimeindx(TCVsig.iteftav)=2;

% initialize order of substructures and allows just a "return" if data empty
trace.data=[];
trace.x=[];
trace.t=[];
trace.dim=[];
trace.dimunits=[];
trace.units=[];

irpintwarn=0;
% find index of signal called upon
if strcmp(data_type_eff(1:1),'\')
  % in case full node name was given
  index=strmatch(data_type_eff(1:end-i_23),TCVsiglocation,'exact');
  if irpintwarn & isempty(index)
    disp('********************')
    disp('trace not yet registered.')
    disp('If standard data, ask andrea.scarabosio@epfl.ch or olivier.sauter@epfl.ch to create a keyqord entry for this data')
%    eval(['!mail -s ''' data_type_eff ' ' num2str(shot) ' ' getenv('USER') ' TCV'' olivier.sauter@epfl.ch < /dev/null'])
    disp('********************')
  elseif isempty(index)
    % temporarily add entry in arrays, so can work below
    index=length(TCVkeywrdall)+1;
    TCVkeywrdall(end+1)={'new'};
    TCVkeywrdcase(end+1)={'simpletdi'};
    TCVsiglocation(end+1)={data_type_eff(1:end-i_23)};
    TCVsigtimeindx(end+1)=0;
  elseif ~strcmp(TCVkeywrdcase{index},'simpletdi')
    msgbox(['Problem in loadTCVdata with data_type_eff = ' data_type_eff ...
          '. Full paths of nodes should only be for case simpletdi'],'in loadTCVdata','error')
    error('in loadTCVdata')
  end
else
  index=strmatch(data_type_eff_noext,TCVkeywrdall,'exact');
  if isempty(index)
    disp(' ')
    disp('********************')
    disp(['no such keyword: ' data_type_eff(1:end-i_23)])
    disp(' ')
    disp('Available keywords:')
    TCVkeywrdall(:)
    disp('********************')
    return
  end
end
if irpintwarn
  disp(['loading' ' ' data_type_eff_noext ' from TCV shot #' num2str(shot)]);
  disp(['case ' TCVkeywrdcase{index}])
end
if i_23==2 & isempty(strmatch(TCVsiglocation(index),liuqe23,'exact')) & strcmp(TCVkeywrdcase{index},'simpletdi')
  disp('********')
  disp('Warning asks for liuqe 2 or 3 of a signal, but not in liuqe23 list in loadTCVdata')
%  eval(['!mail -s ''' data_type_eff ' ' num2str(shot) ' ' getenv('USER') ' TCV: not in liuqe23 list'' olivier.sauter@epfl.ch < /dev/null'])
  disp('********')
end

status=ones(1,100);
zmag=cell(0,0);
nargineff=nargin;

switch TCVkeywrdcase{index}

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case 'simpletdi'

    %  load TCV other data
    mdsopen(shot);
    nodenameeff=[TCVsiglocation{index} liuqe_ext];
    % test if node exists
    error=1;
    ij=findstr(nodenameeff,'[');
    if isempty(ij); ij=length(nodenameeff)+1; end
    if eval(['~mdsdata(''node_exists("\' nodenameeff(1:ij-1) '")'')'])
      disp(['node ' nodenameeff(1:ij-1) ' does not exist for shot = ' num2str(shot)])
      return
% $$$     elseif eval(['mdsdata(''getnci("\' nodenameeff ':foo","length")'')==0'])
% $$$       disp(['no data for node ' nodenameeff ' for shot = ' num2str(shot)])
% $$$       return
    end
    tracetdi=tdi(nodenameeff);
    mdsclose;
    if isempty(tracetdi.data) | isnan(tracetdi.data)
      disp(['node ' nodenameeff ' is empty for shot = ' num2str(shot)])
      return
    end
    trace.data=tracetdi.data;

    if length(tracetdi.dim)==0
      % scalar
      trace.x=[];
      trace.t=[];      
      elseif TCVsigtimeindx(index)>0 | length(tracetdi.dim)==1
      trace.t=tracetdi.dim{max(1,TCVsigtimeindx(index))};
      ix=3-TCVsigtimeindx(index); % works only for 2D arrays
    else
      % if dim=41 assumes is x (usual for TCV)
      if length(tracetdi.dim{1})==41
        ix=1;
        trace.t=tracetdi.dim{2};
      elseif length(tracetdi.dim{2})==41
        ix=2;
        trace.t=tracetdi.dim{1};
      else
        % time array unknow, assumes time array with values having most number of digits
        ab1=num2str(tracetdi.dim{1}-fix(tracetdi.dim{1}));
        ab2=num2str(tracetdi.dim{2}-fix(tracetdi.dim{2}));
        if size(ab1,2)<size(ab2,2); 
          ix=1; 
          trace.t=tracetdi.dim{2};
        else
          ix=2;
          trace.t=tracetdi.dim{1};
        end
        disp(['assumes time array has more digits, so x from dim{' num2str(ix) '}'])
      end
    end
    if length(tracetdi.dim)==2
      trace.x=tracetdi.dim{ix};
      % make sure data is of (x,t)
      if ix==2
        % transpose data
        trace.data=trace.data';
      end
    elseif length(tracetdi.dim)>2
      msgbox(['data more than 2D (data_type_eff=' data_type_eff ...
            '), check how to save it, contact andrea.scarabosio@epfl.ch or olivier.sauter@epfl.ch'],...
          'in simpletdi','warn')
      warning('in simpletdi of loadTCVdata')
    else
      if max(1,TCVsigtimeindx(index))~=1
        disp('Problems in loadTCVdata, max(1,TCVsigtimeindx(index)) should be 1')
      end
    end
    if length(tracetdi.dim)==1
      trace.dim=[{trace.t}];
      trace.dimunits=tracetdi.dimunits(1);
    elseif length(tracetdi.dim)==2
      trace.dim=[{trace.x} ; {trace.t}];
      trace.dimunits=[tracetdi.dimunits(ix) ; tracetdi.dimunits((2-ix)+1)];
    else
      trace.dim=tracetdi.dim;
      trace.dimunits=tracetdi.dimunits;
    end
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    
    trace.name=[num2str(shot) ';' nodenameeff];
    error=0;

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'ne','te'}
    % ne or Te from Thomson data on raw z mesh vs (z,t)
    mdsopen(shot);
    if strcmp(TCVkeywrdcase{index},'ne')
      nodenameeff='\results::thomson:ne';
      tracetdi=tdi(nodenameeff);
      tracestd=tdi('\results::thomson:ne:error_bar');
    else
      nodenameeff='\results::thomson:te';
      tracetdi=tdi(nodenameeff);
      tracestd=tdi('\results::thomson:te:error_bar');
      trace_abs=tdi('\results::thomson:fir_thom_rat');
    end
    trace.data=tracetdi.data'; % Thomson data as (t,z)
    trace.std=tracestd.data';
    trace.name=[num2str(shot) ';' nodenameeff];
    % add correct dimensions
    time=mdsdata('\results::thomson:times');
    z=mdsdata('\diagz::thomson_set_up:vertical_pos');
    trace.dim=[{z};{time}];
    trace.dimunits=[{'Z [m]'} ; {'time [s]'}];
    trace.x=z;
    trace.t=time;
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    mdsclose('mdsip');

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'nerho','terho'}
    % ne or Te from Thomson data on rho=sqrt(psi_normalised) mesh: (rho,t)
    mdsopen(shot);
    time=mdsdata('\results::thomson:times');
    if strcmp(TCVkeywrdcase{index},'nerho')
      nodenameeff='\results::thomson:ne';
      tracetdi=tdi(nodenameeff);
      if isempty(tracetdi.data)
        ishot=mdsopen(shot)
        tracetdi=tdi(nodenameeff)
      end
      nodenameeff='\results::thomson:ne; error_bar ; fir_thom_rat; (ne,std)*fir_thom_rat';
      tracestd=tdi('\results::thomson:ne:error_bar');
      if shot>=23801
        tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
        if isempty(tracefirrat.data)
          disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty')
        end
      else
        tracefirrat=tdi('\results::thomson:fir_thom_rat');
        tracefirrat.dim{1}=time;
      end
      tracefirrat_data=NaN*ones(size(tracetdi.dim{1}));
      if ~isempty(tracefirrat.data)
        itim=iround(time,tracefirrat.dim{1});
        tracefirrat_data(itim)=tracefirrat.data;
      end
    else
      nodenameeff='\results::thomson:te';
      tracetdi=tdi(nodenameeff);
      nodenameeff='\results::thomson:te; error_bar';
      tracestd=tdi('\results::thomson:te:error_bar');
    end
    trace.data=tracetdi.data'; % Thomson data as (t,z)
    trace.std=tracestd.data';
    if strcmp(TCVkeywrdcase{index},'nerho')
      trace.firrat=tracefirrat_data;
      trace.data_abs=trace.data*diag(tracefirrat_data);
      trace.std_abs=trace.std*diag(tracefirrat_data);
    end
    trace.name=[num2str(shot) ';' nodenameeff];
    % add correct dimensions
    % construct rho mesh
    psi_max=tdi(['\results::thomson:psi_max' liuqe_ext]);
    psiscatvol=tdi(['\results::thomson:psiscatvol' liuqe_ext]);
    for ir=1:length(psiscatvol.dim{2})
      rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))';
    end
    trace.dim=[{rho};{time}];
    trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
    trace.x=rho;
    trace.t=time;
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    mdsclose;

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'nerhozshift','terhozshift'}
    % ne or Te from Thomson data on rho=sqrt(psi_normalised) mesh: (rho,t)
    % allow for z shift of equil
    if  nargin>=2 & ~isempty(varargin{1})
      zshift=varargin{1};
    else
      zshift=0.;
    end
    mdsopen(shot);
    if strcmp(TCVkeywrdcase{index},'nerhozshift')
      nodenameeff='\results::thomson:ne';
      tracetdi=tdi(nodenameeff);
      tracestd=tdi('\results::thomson:ne:error_bar');
    else
      nodenameeff='\results::thomson:te';
      tracetdi=tdi(nodenameeff);
      tracestd=tdi('\results::thomson:te:error_bar');
    end
    trace.data=tracetdi.data'; % Thomson data as (t,z)
    trace.std=tracestd.data';
    trace.name=[num2str(shot) ';' nodenameeff];
    % add correct dimensions
    time=mdsdata('\results::thomson:times');
    % construct rho mesh
    psi_max=tdi(['\results::thomson:psi_max' liuqe_ext]);
    psiscatvol=tdi(['\results::thomson:psiscatvol' liuqe_ext]);
    if abs(zshift)>1e-5
      % calculate new psiscatvol
      psitdi=tdi('\results::psi');
      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
          disp(' bad time dimension for zshift')
          disp(['it should be 1 or ' num2str(length(psiscatvol.dim{1})) ' or ' num2str(length(psitdi.dim{3}))])
      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
    for ir=1:length(psiscatvol.dim{2})
      rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))';
    end
    trace.dim=[{rho};{time}];
    trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
    trace.x=rho;
    trace.t=time;
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    mdsclose;

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'profnerho','profterho'}
    % vol from psitbx
    mdsopen(shot);
    if strcmp(TCVkeywrdcase{index},'profnerho')
      nodenameeff=['\results::THOMSON.PROFILES.AUTO:ne'];
      avers=tdi('\results::THOMSON.PROFILES.AUTO:ne:version_num');
    end
    if strcmp(TCVkeywrdcase{index},'profterho')
      nodenameeff=['\results::THOMSON.PROFILES.AUTO:te'];
      avers=tdi('\results::THOMSON.PROFILES.AUTO:te:version_num');
    end
    if avers.data==0
     % may be because nodes not yet filled in, so call once a node
     ab=tdi(nodenameeff);
     avers=tdi([nodenameeff ':version_num']);
    end
    if avers.data>0
      tracetdi=tdi(nodenameeff);
      if avers.data < 2.99
        % for earlier version the bug made it to have logically (rho,t)
        if ~isempty(tracetdi.dim)
          trace.x=tracetdi.dim{1};
          trace.t=tracetdi.dim{2};
        else
          trace.x=[];
          trace.t=[];
        end
      else
        trace.data=tracetdi.data'; % error in dimensions for autofits
        if ~isempty(tracetdi.dim)
          disp('assumes dim{2} for x in THOMSON.PROFILES.AUTO')
          trace.x=tracetdi.dim{2};
          trace.t=tracetdi.dim{1};
        else
          trace.x=[];
          trace.t=[];
        end
      end
    else
      tracetdi=avers;
      trace.x=[];
      trace.t=[];
    end
    trace.dim=[{trace.x};{trace.t}];
    trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    trace.name=[num2str(shot) '; Thomson autfits from ' nodenameeff ];
    mdsclose;

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'qrho'}
    % q profile on psi from liuqe
    mdsopen(shot);
    nodenameeff=['\results::q_psi' liuqe_ext];
    tracetdi=tdi(nodenameeff);
    trace.data=tracetdi.data;
    trace.x=sqrt(tracetdi.dim{1}/(length(tracetdi.dim{1})-1));
    trace.t=tracetdi.dim{2};
    trace.dim=[{trace.x};{trace.t}];
    trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    trace.name=[num2str(shot) ';' nodenameeff];
    mdsclose;

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'vol'}
    % vol from psitbx
    mdsopen(shot);
    nodenameeff=['\results::psitbx:vol'];
    tracetdi=tdi(nodenameeff);
    if i_23==2
      ['LIUQE' liuqe_ext(2:2)]
      psi=psitbxtcv(shot,'+0',['LIUQE' liuqe_ext(2:2)]);
      fsd=psitbxp2p(psi,'FS');
      fsg=psitbxfsg(fsd);
      tracetdi.data = fsg.vol.x;
      tracetdi.dim{1}=fsg.vol.grid.x{:}';
      tracetdi.dim{2}=fsg.vol.grid.t';
    end
    trace.data=tracetdi.data;
    if isempty(tracetdi.data)
      trace.x=tracetdi.dim;
      trace.t=tracetdi.dim;
      trace.dim=tracetdi.dim;
    else
      trace.x=tracetdi.dim{1};
      trace.t=tracetdi.dim{2};
      trace.dim=[{trace.x};{trace.t}];
    end
    trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    trace.name=[num2str(shot) ';' nodenameeff liuqe_ext];
    mdsclose;

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'rhovol'}
    % vol from psitbx
    mdsopen(shot);
    nodenameeff=['\results::psitbx:vol'];
    tracetdi=tdi(nodenameeff);
    if i_23==2
      ['LIUQE' liuqe_ext(2:2)]
      psi=psitbxtcv(shot,'+0',['LIUQE' liuqe_ext(2:2)]);
      fsd=psitbxp2p(psi,'FS');
      fsg=psitbxfsg(fsd);
      tracetdi.data = fsg.vol.x;
      tracetdi.dim{1}=fsg.vol.grid.x{:}';
      tracetdi.dim{2}=fsg.vol.grid.t';
    end
    trace.data=tracetdi.data;
    for i=1:size(tracetdi.data,2)
      trace.data(:,i)=sqrt(tracetdi.data(:,i)./tracetdi.data(end,i));
    end
    trace.x=tracetdi.dim{1};
    trace.t=tracetdi.dim{2};
    trace.dim=[{trace.x};{trace.t}];
    trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    trace.name=[num2str(shot) '; sqrt(V/Va) from ' nodenameeff liuqe_ext];
    mdsclose;

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'sxr','sxR'}
    %  load TCV soft x-ray data
    if  nargin>=3 & ~isempty(varargin{1})
      i1=varargin{1}(1);
      i2=varargin{1}(2);
    else
      i1=1;
      i2=20;
    end
    if  nargin>=4 & ~isempty(varargin{2})
      status=varargin{2};
    else
      status=ones(i2-i1+1,1);
    end

    % camera selection: 1-10. each camera has 20 channels
    icamera=[0 1 0 0 0 0 0 0 0 0];  %index of the camera to use
    if ~isempty(find(status(1:20*icamera*ones(10,1)) == 1))
      [fans,vangle,xchord,ychord,aomega,angfact]=xtomo_geometry(1,icamera);
      % calculating intersection of the view lines with magnetic axis
      if strcmp(data_type_eff,'sxR')
        if nargin>=5 & ~isempty(varargin{3})
          zmag=varargin{3};
        else
          zmag=loadTCVdata(shot,['zmag' liuqe_ext]);
        end
        t_1=zmag.t(1);
        t_2=zmag.t(end);
	[xtomo_signal,t]=get_xtomo_data(shot,t_1,t_2,13e-6*16,icamera,angfact);
	data=interp1(zmag.t,zmag.data,t');
	radius.data=VsxrTCVradius(data,xchord,ychord)';
	radius.t=t';
	varargout{1}={radius};
        trace.R=radius.data;
      else
        t_1=0.001;
        t_2=3;
	[xtomo_signal,t]=get_xtomo_data(shot,t_1,t_2,13e-6*16, ...
          icamera,angfact);
      end
    end
    trace.data=xtomo_signal;
    trace.x=[1:size(trace.data,1)]';
    trace.t=t';
    trace.dim=[{trace.x} ; {trace.t}];
    trace.dimunits=[{'channel #'} ; {'time [s]'}];
    trace.name=[num2str(shot) ';' 'get_xtomo_data'];
    error=0;

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case 'ece'
    %  load TCV ECE data
    % Status=1 => Not Read Yet
    mdsopen(shot);
    if ~isempty(find(status == 1))
      if eval(['~mdsdata(''node_exists("\\RESULTS::ECE:rho")'')'])
        disp(['node \RESULTS::ECE:rho does not exist for shot = ' num2str(shot)])
        return
      end
      if eval(['mdsdata(''getnci("\\RESULTS::ECE:rho","length")'')==0'])
        disp(['no data for \RESULTS::ECE:rho for shot = ' num2str(shot)])
        return
      end
      [TE_ECE,TE_ECE_ERR,RHO,R,T,TE_THOM,TE_THOM_ERR,Fcentral,CAL]=ece_te ...
          (shot,[0.1 0.29],10,10);
    end
    a=min(find(R(:,1)>=0));
    b=max(find(R(:,1)>=0));
    trace.data=TE_ECE(a:b,:)';
    trace.t=T(a:b);
    trace.x=[1:size(trace.data,1)]';
    trace.dim=[{trace.x} ; {trace.t}];
    trace.dimunits=[{'channel #'} ; {'time [s]'}];
    trace.R=R(a:b,:)';
    trace.name=[num2str(shot) ';' '\results::ece...'];
    radius.data=trace.R;
    radius.t=trace.t;
    varargout{1}={radius};
    error=0;
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    mdsclose;
    
  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case 'MPX'
    %  load TCV MPX data
    % Status=1 => Not Read Yet
    if isempty(zmag)
      zmag=loadTCVdata(shot,'zmag');
    end
    t_1=zmag.t(1);
    t_2=zmag.t(end);
    if ~isempty(find(status == 1))
      mdsopen(shot);
      keyboard
      signal=get_mds_mio('MPX',[t_1 t_2]);
      mdsclose;
      trace.data=signal.data;
      for i=1:size(signal.dim{2},2)
        trace.t(:,i)=signal.dim{1};
      end
    end
    trace.dim{1}={trace.t};
    trace.dimunits={'time [s]'};
    trace.name=[num2str(shot) ';' 'get_mds_mio(MPX)'];
    [xchord,ychord]=mpx_geometry;
    varargout{1}={VsxrTCVradius(zmag.data,xchord,ychord)};
    error=0;

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case 'IOH'
    %  Ohmic transformer current
    mdsopen(shot);
    nodenameeff=[{'\magnetics::ipol[*,$1]'} {'OH_001'}];
    tracetdi=tdi(nodenameeff{:});
    trace.data=tracetdi.data;
    trace.x=[];
    trace.t=tracetdi.dim{1};
    trace.dim=tracetdi.dim;
    trace.dimunits={'time [s]'};
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    trace.name=[num2str(shot) ';' nodenameeff{1} ',' nodenameeff{2}];
    mdsclose;
    error=0;
    
  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case 'vloop'
    %  Loop voltage
    mdsopen(shot);
    nodenameeff=[{'\magnetics::vloop[*,$1]'} {'001'}];
    tracetdi=tdi(nodenameeff{:});
    trace.data=tracetdi.data;
    trace.x=[];
    trace.t=tracetdi.dim{1};
    trace.dim=tracetdi.dim;
    trace.dimunits={'time [s]'};
    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work
    if any(strcmp(fieldnames(tracetdi),'units'))
      trace.units=tracetdi.units;
    end
    trace.name=[num2str(shot) ';' nodenameeff{1} ',' nodenameeff{2}];
    mdsclose;
    error=0;
    
  otherwise
    % eval(['!mailto_Andrea ''from loadTCVdata, data_type_eff= ' data_type_eff ''''])
    disp(['this data_type_eff' ' ' data_type_eff ' ' 'not yet programmed in loadTCVdata, ask Andrea.Scarabosio@epfl.ch']);

end