function [trace,error,varargout]=loadAUGdata(shot,data_type,varargin)
%
% data_type:
% 'Ip'   =  current
% 'zmag' =  vertical position of the center of the plasma (magnetic axis)
% 'rmag' =  radial position of the center of the plasma 
% 'sxr'  =  soft x-ray emission 
% 'sxR'  =  soft x-ray emission with varargout{1} option (requires varargin{5}!)
% 'SXB'  =  soft x-ray emission from (by default camera J) SXB/J_xx camera (sxb, sxB, etc all work)
%
%     gdat(15133,'MAG/Ipa',1,'AUG') % uses experiment=AUGD shotfiles per default
%     gdat(15133,'MAG/Ipa',1) (sufficient at AUG since 'AUG' per defaut, same with gdat(15133,'ip',1)
%     gdat(15133,'AUGD/MAG/Ipa',1,'AUG') % to specify experiment explicitely like in:
%     gdat(30230,'ECED/RMD/Trad-A',1,'AUG') % 
%
% INPUT:
% shot: shot number
% data_type: type of the required data: 'diag_name/sig_name'
%
% examples:
%           data_type='SXR/B', 'TOT/beta_N', 'SXB/J_053', 'SXB'
%           data_type='POT/ELMa-Han', 'MOD/OddNAmp', 'MOD/EvenNAmp', 'TOT/PNBI_TOT', 'TOT/P_TOT'
%
% Meaning of varargin depends on data_type:
%
% data_type=sxr, sxb or ece, eced:
%                  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
%                  varargin{4}:  time range [t1 t2] (to limit data collected)
%                  varargin{5}:  nth to keep only nth time points of traces
%                  varargin{6}:  camera to use: 'B' (default), 'A', 
%
% OUTPUT:
% trace.data:   data structure
% trace.t:      time of reference 
% trace.x:      space of reference 
% ....          others related to data
% error:        error in loading signal (0=> OK, 1=> error)
%
% Additional Output arguments depending on data_type
%
% data_type=sxR:
%                varargout{1}: intersection of the view lines with magnetic axis
%
% functions needed: SF2ML or mdsplus routines
%
% Example:
%         [ip,error]=loadAUGdata(shot,'ip');
%         [ip,error]=loadAUGdata(shot,'MAG/Ipi');
%         [n2,error]=loadAUGdata(shot,'MOD/EvenNAmp');
%
% List of keywords (can be used in adition to 'DIAG/tracenam'), with comments when not obvious:
%
% Ip : 
% b0 : 
% zmag : 
% rmag : 
% rgeo : 
% zgeo : 
% q0 : 
% q95 : 
% kappa : 
% delta : 
% deltatop : 
% deltabot : 
% neint : 
% neterho : 
% cxrs : CXRS structure from CEZ with vrot, Ti, ...
% cxrs_rho : same as cxrs but project on rho as well (using 'equil' information)
% equil : equilibrium structure, rhopol, rhotor, rhovol, etc from EQI
% equil_fpp : as equil but from FPP
% equil_eqm : as equil but from EQM
% equil_eqr : as equil but from EQR
% sxr : from SXR/A or B (from old stuff, not sure still OK)
% sxR : from SXR/A or B adding R of chords (from old stuff, not sure still OK)
% sxb : 'SXB/J' chords
% ece : 
% eced : 
% Halpha : 
% pgyro : for each gyrotrons, power, freq, etc (ask for more)
% powers : power traces for each sources
%


varargout={cell(1,1)};
error=1;
shotfile_exp = 'AUGD';

% To allow multiple ways of writing a specific keyword, use data_type_eff within this routine
data_type_eff=data_type;
if size(data_type_eff,1)==1
  i=findstr('/',data_type);
  if length(i)==1
    % assumes given a la 'MAG/Ipi'
    data_type_eff=[{data_type(1:i(1)-1)} ; {data_type(i(1)+1:end)}];
  elseif length(i)==2
    % assumes given a la 'AUGD/MAG/Ipi' or 'ECED/RMP/Trad-A'
    data_type_eff=[{data_type(i(1)+1:i(2)-1)} ; {data_type(i(2)+1:end)}];
    shotfile_exp = data_type(1:i(1)-1);
  elseif length(i)>2
    disp(['more / than expected in tracename: length(i)= ' num2str(length(i))])
    data_type
  end
end

i_efitm=0;
i_ext=length(data_type_eff)+1;
name_ext='';
if size(data_type_eff,1)==1
  data_type_eff_noext=data_type_eff(1:i_ext-1);
  if ~isempty(strmatch(data_type_eff_noext,[{'ip'} {'i_p'} {'IP'} {'iP'} {'xip'}],'exact'))
    data_type_eff_noext='Ip';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'B0'}],'exact'))
    data_type_eff_noext='b0';
  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,[{'Terho'}],'exact'))
    data_type_eff_noext='terho';
  end
  if ~isempty(strmatch(lower(data_type_eff_noext),[{'neterho'}],'exact'))
    data_type_eff_noext='neterho';
  end
  if ~isempty(strmatch(lower(data_type_eff_noext),[{'cxrs'} {'vrot'} {'ti'}],'exact'))
    data_type_eff_noext='cxrs'; % load full CEZ structure
  end
  if ~isempty(strmatch(lower(data_type_eff_noext),[{'cxrs_rhos'} {'cxrs_rho'} {'cxrsrho'} {'cxrsrhos'}],'exact'))
    data_type_eff_noext='cxrs_rho'; % load full CEZ structure
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'SXR'}],'exact'))
    data_type_eff_noext='sxr';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'SXB'} {'sxb'} {'Sxb'} {'sXb'} {'sxB'} {'SXb'}],'exact'))
    data_type_eff_noext='sxb';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'ECE'}],'exact'))
    data_type_eff_noext='ece';
  end
  if ~isempty(strmatch(upper(data_type_eff_noext),[{'ECED'}],'exact'))
    data_type_eff_noext='eced';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'VOL'} {'volume'}],'exact'))
    data_type_eff_noext='vol';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'q_0'} {'Q0'}],'exact'))
    data_type_eff_noext='q0';
  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,[{'qrho_FPP'}],'exact'))
    data_type_eff_noext='qrho_fpp';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'equil_FPP'}],'exact'))
    data_type_eff_noext='equil_fpp';
  end
  if ~isempty(strmatch(lower(data_type_eff_noext),[{'equil_eqm'}],'exact'))
    data_type_eff_noext='equil_eqm';
  end
  if ~isempty(strmatch(lower(data_type_eff_noext),[{'equil_eqr'}],'exact'))
    data_type_eff_noext='equil_eqr';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'Rmag'}],'exact'))
    data_type_eff_noext='rmag';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'Rgeo'}],'exact'))
    data_type_eff_noext='rgeo';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'Zmag'}],'exact'))
    data_type_eff_noext='zmag';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'Zgeo'}],'exact'))
    data_type_eff_noext='zgeo';
  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,[{'Ha'} {'ha'} {'Halpha'} {'halpha'}],'exact'))
    data_type_eff_noext='Halpha';
  end
  if ~isempty(strmatch(lower(data_type_eff_noext),[{'pgyro'} {'pec'} {'pech'} {'pecrh'} {'p_ec'} {'p_gyro'}],'exact'))
    data_type_eff_noext='pgyro';
  end
  if ~isempty(strmatch(lower(data_type_eff_noext),[{'powers'} {'ptot'} {'ptots'}],'exact'))
    data_type_eff_noext='powers';
  end
else
  i_ext=length(data_type_eff{2})+1;
  name_ext='';
  data_type_eff_noext=data_type_eff{2}(1:i_ext-1);
end

% all keywords and corresponding case to run below
AUGkeywrdall=[{'Ip'} {'b0'} {'zmag'} {'rmag'}  {'rgeo'} {'zgeo'} {'rcont'} {'zcont'} {'vol'} {'qrho'} {'qrho_fpp'} {'q0'} {'q95'} {'kappa'} ...
	      {'delta'} {'deltatop'} {'deltabot'} {'neint'} {'ne'} {'te'}  ...
	      {'nerho'} {'neterho'} {'terho'} {'cxrs'} {'cxrs_rho'} {'equil'} {'equil_fpp'} {'equil_eqm'} ...
	      {'equil_eqr'} {'sxr'} {'sxR'} {'sxb'} {'ece'} {'eced'} {'Halpha'} {'pgyro'} {'powers'}];
AUGsig.iip=strmatch('Ip',AUGkeywrdall,'exact');
AUGsig.ib0=strmatch('b0',AUGkeywrdall,'exact');
AUGsig.izmag=strmatch('zmag',AUGkeywrdall,'exact');
AUGsig.irmag=strmatch('rmag',AUGkeywrdall,'exact');
AUGsig.irgeo=strmatch('rgeo',AUGkeywrdall,'exact');
AUGsig.izgeo=strmatch('zgeo',AUGkeywrdall,'exact');
AUGsig.ircont=strmatch('rcont',AUGkeywrdall,'exact');
AUGsig.izcont=strmatch('zcont',AUGkeywrdall,'exact');
AUGsig.ivol=strmatch('vol',AUGkeywrdall,'exact');
AUGsig.iqrho=strmatch('qrho',AUGkeywrdall,'exact');
AUGsig.iqrho_fpp=strmatch('qrho_fpp',AUGkeywrdall,'exact');
AUGsig.iequil=strmatch('equil',AUGkeywrdall,'exact');
AUGsig.iequil_fpp=strmatch('equil_fpp',AUGkeywrdall,'exact');
AUGsig.iequil_eqm=strmatch('equil_eqm',AUGkeywrdall,'exact');
AUGsig.iequil_eqr=strmatch('equil_eqr',AUGkeywrdall,'exact');
AUGsig.iq0=strmatch('q0',AUGkeywrdall,'exact');
AUGsig.iq95=strmatch('q95',AUGkeywrdall,'exact');
AUGsig.ikappa=strmatch('kappa',AUGkeywrdall,'exact');
AUGsig.idelta=strmatch('delta',AUGkeywrdall,'exact');
AUGsig.ideltatop=strmatch('deltatop',AUGkeywrdall,'exact');
AUGsig.ideltabot=strmatch('deltabot',AUGkeywrdall,'exact');
AUGsig.ineint=strmatch('neint',AUGkeywrdall,'exact');
AUGsig.ine=strmatch('ne',AUGkeywrdall,'exact');
AUGsig.ite=strmatch('te',AUGkeywrdall,'exact');
AUGsig.inerho=strmatch('nerho',AUGkeywrdall,'exact');
AUGsig.iterho=strmatch('terho',AUGkeywrdall,'exact');
AUGsig.ineterho=strmatch('neterho',AUGkeywrdall,'exact');
AUGsig.icxrs=strmatch('cxrs',AUGkeywrdall,'exact');
AUGsig.icxrs_rho=strmatch('cxrs_rho',AUGkeywrdall,'exact');
AUGsig.isxr=strmatch('sxr',AUGkeywrdall,'exact');
AUGsig.isxR=strmatch('sxR',AUGkeywrdall,'exact');
AUGsig.isxb=strmatch('sxb',AUGkeywrdall,'exact');
AUGsig.iece=strmatch('ece',AUGkeywrdall,'exact');
AUGsig.ieced=strmatch('eced',AUGkeywrdall,'exact');
AUGsig.iHalpha=strmatch('Halpha',AUGkeywrdall,'exact');
AUGsig.ipgyro=strmatch('pgyro',AUGkeywrdall,'exact');
AUGsig.ipowers=strmatch('powers',AUGkeywrdall,'exact');

% For each keyword, specify which case to use. As most common is 'simplereaddata', fill in with this and change
% only indices needed. Usually use name of case same as keyword name
AUGkeywrdcase=cell(size(AUGkeywrdall));
AUGkeywrdcase(:)={'simplereaddata'};
AUGkeywrdcase(AUGsig.iqrho)=AUGkeywrdall(AUGsig.iqrho); % special as efit q on psi
AUGkeywrdcase(AUGsig.iqrho_fpp)=AUGkeywrdall(AUGsig.iqrho_fpp); % special as efit q on psi
AUGkeywrdcase(AUGsig.iequil)=AUGkeywrdall(AUGsig.iequil); % special as efit q on psi
AUGkeywrdcase(AUGsig.iequil_fpp)=AUGkeywrdall(AUGsig.iequil_fpp); % special as efit q on psi
AUGkeywrdcase(AUGsig.iequil_eqm)=AUGkeywrdall(AUGsig.iequil_eqm); % special as efit q on psi
AUGkeywrdcase(AUGsig.iequil_eqr)=AUGkeywrdall(AUGsig.iequil_eqr); % special as efit q on psi
%AUGkeywrdcase(AUGsig.idelta)=AUGkeywrdall(AUGsig.idelta); % special as average of triu and tril
AUGkeywrdcase(AUGsig.ine)=AUGkeywrdall(AUGsig.ine); % special as adds error bars
AUGkeywrdcase(AUGsig.ite)=AUGkeywrdall(AUGsig.ite); % idem
AUGkeywrdcase(AUGsig.inerho)=AUGkeywrdall(AUGsig.inerho); % idem
AUGkeywrdcase(AUGsig.ineterho)=AUGkeywrdall(AUGsig.ineterho); % idem
AUGkeywrdcase(AUGsig.iterho)=AUGkeywrdall(AUGsig.iterho); % idem
AUGkeywrdcase(AUGsig.isxr)=AUGkeywrdall(AUGsig.isxr);
AUGkeywrdcase(AUGsig.isxR)=AUGkeywrdall(AUGsig.isxR);
AUGkeywrdcase(AUGsig.isxb)=AUGkeywrdall(AUGsig.isxb);
%AUGkeywrdcase(AUGsig.iece)=AUGkeywrdall(AUGsig.iece);
AUGkeywrdcase(AUGsig.icxrs)=AUGkeywrdall(AUGsig.icxrs);
AUGkeywrdcase(AUGsig.icxrs_rho)=AUGkeywrdall(AUGsig.icxrs_rho);
AUGkeywrdcase(AUGsig.ipgyro)=AUGkeywrdall(AUGsig.ipgyro); % idem
AUGkeywrdcase(AUGsig.ipowers)=AUGkeywrdall(AUGsig.ipowers); % idem

% 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
AUGsigtimeindx=ones(size(AUGkeywrdall));

% For the 'simplereaddata' cases, we need the full node in case gdat was called with full location directly
% for the other cases, leave this location empty
AUGsiglocation=cell(2,size(AUGkeywrdall,2));
AUGsiglocation(:)={''};
AUGsiglocation(:,AUGsig.iip)={'MAG'; 'Ipa'};
AUGsiglocation(:,AUGsig.ib0)={'FPC'; 'BTF'}; % at 1.65m (?)
AUGsiglocation(:,AUGsig.izmag)={'FPG'; 'Zmag'};
AUGsiglocation(:,AUGsig.irmag)={'FPG'; 'Rmag'};
AUGsiglocation(:,AUGsig.irgeo)={'FPG'; 'Rgeo'};
AUGsiglocation(:,AUGsig.izgeo)={'FPG'; 'Zgeo'};
AUGsiglocation(:,AUGsig.ircont)={'' ; ''}; AUGsigtimeindx(AUGsig.ircont)=2;
AUGsiglocation(:,AUGsig.izcont)={'' ; ''}; AUGsigtimeindx(AUGsig.izcont)=2;
AUGsiglocation(:,AUGsig.ivol)={''; ''};
AUGsiglocation(:,AUGsig.iq0)={'FPG'; 'q0'};
AUGsiglocation(:,AUGsig.iq95)={'FPG'; 'q95'};
AUGsiglocation(:,AUGsig.ikappa)={''; ''};
AUGsiglocation(:,AUGsig.ideltatop)={''; ''};
AUGsiglocation(:,AUGsig.ideltabot)={''; ''};
AUGsiglocation(:,AUGsig.ineint)={'DCN'; 'H-1'};
AUGsiglocation(:,AUGsig.iece)={'CEC'; 'Trad-A'};
AUGsiglocation(:,AUGsig.ieced)={'RMD'; 'Trad-A'}; % ECED
AUGsiglocation(:,AUGsig.iHalpha)={'POT'; 'ELMa-Han'};

% For the 'simplereaddata' cases, we need the full node in case gdat was called with full location directly
% for the other cases, leave this location empty
AUGexplocation=cell(size(AUGkeywrdall,2),1);
AUGexplocation(:)={'AUGD'};
% cases with a "private" shotfile:
AUGexplocation(AUGsig.ieced)={'ECED'};

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

% find index of signal called upon
if size(data_type_eff,1)==2
  % in case node name was given in 2 parts directly (as old way)
  ii1=strmatch(data_type_eff(1),AUGsiglocation(1,:),'exact');
  iiindex=strmatch(data_type_eff_noext,AUGsiglocation(2,ii1),'exact');
  if ~isempty(iiindex)
    index=ii1(iiindex);
  else
    index=[];
  end
  if 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{1} ' ' data_type_eff{2} ' ' num2str(shot) ' ' ...
%          getenv('USER') ' AUG'' olivier.sauter@epfl.ch < /dev/null'])
    disp('********************')
    % temporarily add entry in arrays, so can work below
    index=length(AUGkeywrdall)+1;
    AUGkeywrdall(end+1)={'new'};
    AUGkeywrdcase(end+1)={'simplereaddata'};
    AUGsiglocation(1:2,end+1)=[data_type_eff(1) ; {data_type_eff_noext}];
    AUGexplocation{end+1}=shotfile_exp;
    AUGsigtimeindx(end+1)=0;
  elseif ~strcmp(AUGkeywrdcase{index},'simplereaddata')
    msgbox(['Problem in loadAUGdata with data_type_eff = ' char(data_type_eff(end)) ...
          '. Full paths of nodes should only be for case simplereaddata'],'in loadAUGdata','error')
    error('in loadAUGdata')
  end
else
  index=strmatch(data_type_eff_noext,AUGkeywrdall,'exact');
  if isempty(index)
    disp(' ')
    disp('********************')
    if iscell(data_type_eff)
      disp(['no such keyword: ' data_type_eff{1} '/' data_type_eff{2}])
    else
      disp(['no such keyword: ' data_type_eff])
    end
    disp(' ')
    disp('Available keywords:')
    AUGkeywrdall(:)
    disp('********************')
    return
  end
end
disp(' ')
if iscell(data_type_eff)
  disp(['loading' ' ' data_type_eff{1} '/' data_type_eff{2} ' from AUG shot #' num2str(shot)]); 
else
  disp(['loading' ' ' data_type_eff ' from AUG shot #' num2str(shot)]); 
end
disp(['case ' AUGkeywrdcase{index}])
disp(' ')
switch AUGkeywrdcase{index}

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case 'simplereaddata'

    ppftype=AUGsiglocation{1,index};
    shotfile_exp_eff = AUGexplocation{index};
    if i_efitm; 
      tracename=['eftm' AUGsiglocation{2,index}(5:end) name_ext];
    else
      tracename=[AUGsiglocation{2,index} name_ext];
    end
    ij=find(tracename~='''');
    tracename=tracename(ij)
    [a,e]=rdaAUG_eff(shot,ppftype,tracename,shotfile_exp_eff);
%    switch tracename
%  special cases if traces do not exist for some shot or other
%    end
    
    trace=a;
    clear error
    error=e;
    if isempty(trace.data)
      trace.dimunits=[];
    elseif length(size(trace.data))==1 | (length(size(trace.data))==2 & size(trace.data,2)==1)
      trace.dim=[{trace.t}];
      trace.dimunits={'time [s]'};
    elseif length(size(trace.data))==2
      trace.dim=[{trace.x} ; {trace.t}];
      trace.dimunits=[{'R [m] or rho=sqrt(psi_norm)'} ; {'time [s]'}];
    else
      disp('how to deal with 3D arrays?')
      trace.dim=[{trace.x} ; {trace.t} ; {d}];
      trace.dimunits=[{'R [m] or rho=sqrt(psi_norm)'} ; {'time [s]'} ; {'d'}];
      trace.d=d;
    end
    trace.name=[num2str(shot) '/' ppftype '/' tracename];
    
  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'sxr','sxR'}
    %  LOAD MULTI CHANNEL DATA    
    %  load AUG soft x-ray data 

    shotfile_exp_eff = AUGexplocation{index};

    if  nargin>=3 & ~isempty(varargin{1})
      starti=varargin{1}(1);
      endi=varargin{1}(2);
    else
      starti=1;
      endi=30;
    end
    if  nargin>=4 & ~isempty(varargin{2})
      status=varargin{2};
    else
      status=ones(endi-starti+1,1);
    end
    if  nargin>=6 & ~isempty(varargin{4})
      timerange=varargin{4};
    else
      timerange=[1 7];
    end
    if  nargin>=7 & ~isempty(varargin{5})
      nth=varargin{5};
    else
      nth=13;
    end
    if  nargin>=8 & ~isempty(varargin{6})
      tracename=varargin{6};
    else
      tracename='B';
    end
    trace.t=[];
    trace.x=[];
    ppftype='SXR';
    [a,e]=rdaAUG_eff(shot,ppftype,tracename,shotfile_exp_eff,timerange);
    trace=a;
    trace.dim=[{[starti:endi]'} ; {trace.t}];
    trace.x=trace.dim{1};
    trace.dimunits=[{'channels'} ; {'time [s]'}];
    trace.units='W/m^2';
    trace.name=[num2str(shot) '/' ppftype '/' tracename];
    % keep only nth points
    trace.t=trace.t(1:nth:end);
    trace.data=trace.data(:,1:nth:end);
    trace.dim{2}=trace.t;
    % calculating intersection of the view lines with magnetics axis 
    if strcmp(data_type_eff_noext,'sxR')
      if nargin>=5 & ~isempty(varargin{3})
        zmag=varargin{3};
      else
        zmag=loadAUGdata(shot,'zmag');
      end
      zmageff=interp1(zmag.t,zmag.data,trace.t);
      if strcmp(tracename,'B')
        [R_B, Z_B, ang_B,Rsxr]=sxrbgeometry(zmageff);
      elseif strcmp(tracename,'A')
        [R_A, Z_A, ang_A,Rsxr]=sxrageometry(zmageff);
      else
        disp(['sxr camera: ' tracename ' not set yet for calculating R projection'])
        return
      end
      radius.data=Rsxr;
      radius.t=trace.t;
      varargout{1}={radius};
      trace.R=radius;
    end

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'sxb'}
    %  LOAD MULTI CHANNEL DATA SXB/J_0xx (or other than J camera if specified in varargin{8})
    %  load AUG soft x-ray data 

    shotfile_exp_eff = AUGexplocation{index};

    if  nargin>=3 & ~isempty(varargin{1})
      % chords to be loaded
      starti=varargin{1}(1);
      endi=varargin{1}(2);
    else
      starti=52;
      endi=54;
    end
    if  nargin>=4 & ~isempty(varargin{2})
      % chords already loaded, 1=not loaded (to load), 0=loaded
      status=varargin{2};
    else
      status=ones(endi-starti+1,1);
    end
    if  nargin>=6 & ~isempty(varargin{4})
      timerange=varargin{4};
    else
      timerange=[0 10];
    end
    if  nargin>=7 & ~isempty(varargin{5})
      nth=varargin{5};
    else
      nth=13;
    end
    if  nargin>=8 & ~isempty(varargin{6})
      tracename=varargin{6};
    else
      tracename='J';
    end
    trace.t=[];
    trace.x=[];
    ppftype='SXB';
    iok=0;
    for ichord=starti:endi
      tracename_eff = [tracename '_' num2str(ichord,'%.3d')];
      [a,e]=rdaAUG_eff(shot,ppftype,tracename_eff,shotfile_exp_eff,timerange);
      if isempty(a) || e~=0
	trace_all = struct([]);
      else
	if iok==0
	  trace_all = a;
	  trace_all = rmfield(trace_all,[{'value'},{'data'}]);
	  iok = iok+1;
	else
	  iok = iok+1;
	end
	trace_all.value(ichord,:) = a.value;
	trace_all.data(ichord,:) = a.data;
      end
    end
    if ~isempty(trace_all)
      trace_all.dim=[{[starti:endi]'} ; {trace.t}];
      trace = trace_all;
      trace.x=trace.dim{1};
      trace.dimunits=[{'channels'} ; {'time [s]'}];
      trace.units='W/m^2';
      trace.name=[num2str(shot) '/' ppftype '/' tracename];
      % keep only nth points
      trace.t=trace.t(1:nth:end);
      trace.data=trace.data(:,1:nth:end);
      trace.dim{2}=trace.t;
      trace.value=trace.value(:,1:nth:end);
      trace.time_aug.value=trace.time_aug.value(1:nth:end);
    else
      trace.data = [];
      trace.dim = [];
      trace.dimunits = [];
      trace.x = [];
      trace.t = [];
      trace.units = [];
      trace.name=[num2str(shot) '/' ppftype '/' tracename];
    end
    
  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'cxrs', 'cxrs_rho'}
    %  LOAD CEZ data
    %  

    shotfile_exp_eff = AUGexplocation{index};

    trace.t=[];
    trace.x=[];
    ppftype='CEZ';
    [a,e]=rdaAUG_eff(shot,ppftype,'vrot',shotfile_exp_eff);
    if isempty(a) || e~=0
      cxrs = struct([]);
    else
      cxrs.vrot = a;
      rmfield(cxrs.vrot,'x'); rmfield(cxrs.vrot,'t'); rmfield(cxrs.vrot,'data');
      cxrs.x = a.x;
      cxrs.data = a.data;
      cxrs.t = a.t;
      cxrs.x = a.x;
    end
    if ~isempty(cxrs)
      cxrs.dim=[{cxrs.x} ; {cxrs.t}];
      cxrs.dimunits=[{'chord'} ; {'time [s]'}];
      cxrs.name=[num2str(shot) '/' ppftype '/vrot;Ti;Ti_c'];
      [aerr,e]=rdaAUG_eff(shot,ppftype,'err_vrot',shotfile_exp_eff);
      cxrs.vrot.error = aerr.value;
      [r_time]=sf2ab(ppftype,shot,'R_time','-exp',shotfile_exp_eff);
      cxrs.r_time = r_time.value{1};
      [z_time]=sf2ab(ppftype,shot,'z_time','-exp',shotfile_exp_eff);
      cxrs.z_time = z_time.value{1};
      [a,e]=rdaAUG_eff(shot,ppftype,'Ti',shotfile_exp_eff);
      [aerr,e]=rdaAUG_eff(shot,ppftype,'err_Ti',shotfile_exp_eff);
      cxrs.ti = a;
      cxrs.ti.error = aerr.value;
      [a,e]=rdaAUG_eff(shot,ppftype,'Ti_c',shotfile_exp_eff);
      [aerr,e]=rdaAUG_eff(shot,ppftype,'err_Ti_c',shotfile_exp_eff);
      cxrs.ti_c = a;
      cxrs.ti_c.error = aerr.value;
      %
      if strcmp(AUGkeywrdcase{index},'cxrs_rho')
	equil=gdat(shot,'equil',0);
	inb_chord_cxrs=size(cxrs.r_time,1);
	inb_time_cxrs=size(cxrs.r_time,2);
	psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
	rhopsinorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
	rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
	rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
	% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
	time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-1)];
	iok=find(cxrs.r_time(:,1)>0);
	for itequil=1:length(time_equil)-1
	  rr=equil.Rmesh(:,itequil);
	  zz=equil.Zmesh(:,itequil);
	  psirz_in = equil.psi2D(:,:,itequil);
	  it_cxrs_inequil = find(cxrs.t>=time_equil(itequil) & cxrs.t<=time_equil(itequil+1));
	  if ~isempty(it_cxrs_inequil)
	    rout=cxrs.r_time(iok,it_cxrs_inequil);
	    zout=cxrs.z_time(iok,it_cxrs_inequil);
	    psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
	    psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil));
	    rhopsinorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
	    for it_cx=1:length(it_cxrs_inequil)
	      rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
	      rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
	    end
	  end
	end
	cxrs.psi_on_rztime = psi_out;
	cxrs.rhopsinorm_on_rztime = rhopsinorm_out;
	cxrs.rhotornorm_on_rztime = rhotornorm_out;
	cxrs.rhovolnorm_on_rztime = rhovolnorm_out;
      end
	
      trace = cxrs;
    else
      trace.data = [];
      trace.dim = [];
      trace.dimunits = [];
      trace.x = [];
      trace.t = [];
      trace.units = [];
      trace.name=[num2str(shot) '/' ppftype '/' tracename];
    end
    
  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'neterho'} % 'nerho', 'terho', 
    %  LOAD VTA data
    %  Vertical Thomson core and edge

    shotfile_exp_eff = AUGexplocation{index};

    trace.t=[];
    trace.x=[];
    ppftype='VTA';
    % if strcmp(AUGkeywrdcase{index},'terho')    
    [a,e]=rdaAUG_eff(shot,ppftype,'Te_c',shotfile_exp_eff);
    if isempty(a) || e~=0
      thomson = struct([]);
    else
      thomson.te_core = a;
      rmfield(thomson.te_core,'x'); rmfield(thomson.te_core,'t'); rmfield(thomson.te_core,'data');
      thomson.x = a.x;
      thomson.data = a.data;
      thomson.t = a.t;
      thomson.time_core = a.t;
      thomson.x = a.x;
    end
    if ~isempty(thomson)
      thomson.dim=[{thomson.x} ; {thomson.t}];
      thomson.dimunits=[{'points'} ; {'time [s]'}];
      thomson.name=[num2str(shot) '/' ppftype '/te_c;te_e;ne_c;ne_e'];
      [alow,e]=rdaAUG_eff(shot,ppftype,'Telow_c',shotfile_exp_eff);
      [aup,e]=rdaAUG_eff(shot,ppftype,'Teupp_c',shotfile_exp_eff);
      thomson.te_core.error = max(aup.value-thomson.te_core.value,thomson.te_core.value-alow.value);
      %
      [a,e]=rdaAUG_eff(shot,ppftype,'Ne_c',shotfile_exp_eff);
      thomson.ne_core = a;
      rmfield(thomson.ne_core,'x'); rmfield(thomson.ne_core,'t'); rmfield(thomson.ne_core,'data');
      [alow,e]=rdaAUG_eff(shot,ppftype,'Nelow_c',shotfile_exp_eff);
      [aup,e]=rdaAUG_eff(shot,ppftype,'Neupp_c',shotfile_exp_eff);
      thomson.ne_core.error = max(aup.value-thomson.ne_core.value,thomson.ne_core.value-alow.value);
      %
      [r_time]=rdaAUG_eff(shot,ppftype,'R_core',shotfile_exp_eff);
      thomson.r_core_time = ones(size(thomson.te_core.value,1),1)*r_time.value;
      [z_time]=rdaAUG_eff(shot,ppftype,'Z_core',shotfile_exp_eff);
      thomson.z_core_time = z_time.value' * ones(1,size(thomson.te_core.value,2));
      %
      [a,e]=rdaAUG_eff(shot,ppftype,'Te_e',shotfile_exp_eff);
      thomson.te_edge = a;
      thomson.time_edge = a.t;
      rmfield(thomson.te_edge,'x'); rmfield(thomson.te_edge,'t'); rmfield(thomson.te_edge,'data');
      [alow,e]=rdaAUG_eff(shot,ppftype,'Telow_e',shotfile_exp_eff);
      [aup,e]=rdaAUG_eff(shot,ppftype,'Teupp_e',shotfile_exp_eff);
      thomson.te_edge.error = max(aup.value-thomson.te_edge.value,thomson.te_edge.value-alow.value);
      %
      [a,e]=rdaAUG_eff(shot,ppftype,'Ne_e',shotfile_exp_eff);
      thomson.ne_edge = a;
      rmfield(thomson.ne_edge,'x'); rmfield(thomson.ne_edge,'t'); rmfield(thomson.ne_edge,'data');
      [alow,e]=rdaAUG_eff(shot,ppftype,'Nelow_e',shotfile_exp_eff);
      [aup,e]=rdaAUG_eff(shot,ppftype,'Neupp_e',shotfile_exp_eff);
      thomson.ne_edge.error = max(aup.value-thomson.ne_edge.value,thomson.ne_edge.value-alow.value);
      %
      [r_time]=rdaAUG_eff(shot,ppftype,'R_edge',shotfile_exp_eff);
      thomson.r_edge_time = ones(size(thomson.te_edge.value,1),1)*r_time.value;
      [z_time]=rdaAUG_eff(shot,ppftype,'Z_edge',shotfile_exp_eff);
      thomson.z_edge_time = z_time.value' * ones(1,size(thomson.te_edge.value,2));
      %
      if strcmp(AUGkeywrdcase{index},'neterho')
	equil=gdat(shot,'equil',0);
	% core
	inb_chord_thomson_core=size(thomson.r_core_time,1);
	inb_time_thomson_core=size(thomson.r_core_time,2);
	psi_out_core = NaN*ones(inb_chord_thomson_core,inb_time_thomson_core);
	rhopsinorm_out_core = NaN*ones(inb_chord_thomson_core,inb_time_thomson_core);
	rhotornorm_out_core = NaN*ones(inb_chord_thomson_core,inb_time_thomson_core);
	rhovolnorm_out_core = NaN*ones(inb_chord_thomson_core,inb_time_thomson_core);
	% edge
	inb_chord_thomson_edge=size(thomson.r_edge_time,1);
	inb_time_thomson_edge=size(thomson.r_edge_time,2);
	psi_out_edge = NaN*ones(inb_chord_thomson_edge,inb_time_thomson_edge);
	rhopsinorm_out_edge = NaN*ones(inb_chord_thomson_edge,inb_time_thomson_edge);
	rhotornorm_out_edge = NaN*ones(inb_chord_thomson_edge,inb_time_thomson_edge);
	rhovolnorm_out_edge = NaN*ones(inb_chord_thomson_edge,inb_time_thomson_edge);
	% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
	time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-1)];
	for itequil=1:length(time_equil)-1
	  rr=equil.Rmesh(:,itequil);
	  zz=equil.Zmesh(:,itequil);
	  psirz_in = equil.psi2D(:,:,itequil);
	  it_thomson_core_inequil = find(thomson.time_core>=time_equil(itequil) & thomson.time_core<=time_equil(itequil+1));
	  if ~isempty(it_thomson_core_inequil)
	    rout_core=thomson.r_core_time(:,it_thomson_core_inequil);
	    zout_core=thomson.z_core_time(:,it_thomson_core_inequil);
	    psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core);
	    psi_out_core(:,it_thomson_core_inequil) = reshape(psi_at_routzout,inb_chord_thomson_core,length(it_thomson_core_inequil));
	    rhopsinorm_out_core(:,it_thomson_core_inequil) = sqrt((psi_out_core(:,it_thomson_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
	    for it_cx=1:length(it_thomson_core_inequil)
	      rhotornorm_out_core(:,it_thomson_core_inequil(it_cx)) = ...
		  interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out_core(:,it_thomson_core_inequil(it_cx)),-3,[2 2],[0 1]);
	      rhovolnorm_out_core(:,it_thomson_core_inequil(it_cx)) = ...
		  interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out_core(:,it_thomson_core_inequil(it_cx)),-3,[2 2],[0 1]);
	    end
	  end
	  % edge
	  it_thomson_edge_inequil = find(thomson.time_edge>=time_equil(itequil) & thomson.time_edge<=time_equil(itequil+1));
	  if ~isempty(it_thomson_edge_inequil)
	    rout_edge=thomson.r_edge_time(:,it_thomson_edge_inequil);
	    zout_edge=thomson.z_edge_time(:,it_thomson_edge_inequil);
	    psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_edge,zout_edge);
	    psi_out_edge(:,it_thomson_edge_inequil) = reshape(psi_at_routzout,inb_chord_thomson_edge,length(it_thomson_edge_inequil));
	    rhopsinorm_out_edge(:,it_thomson_edge_inequil) = sqrt((psi_out_edge(:,it_thomson_edge_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
	    for it_cx=1:length(it_thomson_edge_inequil)
	      rhotornorm_out_edge(:,it_thomson_edge_inequil(it_cx)) = ...
		  interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out_edge(:,it_thomson_edge_inequil(it_cx)),-3,[2 2],[0 1]);
	      rhovolnorm_out_edge(:,it_thomson_edge_inequil(it_cx)) = ...
		  interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out_edge(:,it_thomson_edge_inequil(it_cx)),-3,[2 2],[0 1]);
	    end
	  end
	end
	thomson.core_psi_on_rztime = psi_out_core;
	thomson.core_rhopsinorm_on_rztime = rhopsinorm_out_core;
	thomson.core_rhotornorm_on_rztime = rhotornorm_out_core;
	thomson.core_rhovolnorm_on_rztime = rhovolnorm_out_core;
	thomson.edge_psi_on_rztime = psi_out_edge;
	thomson.edge_rhopsinorm_on_rztime = rhopsinorm_out_edge;
	thomson.edge_rhotornorm_on_rztime = rhotornorm_out_edge;
	thomson.edge_rhovolnorm_on_rztime = rhovolnorm_out_edge;
      end
	
      trace = thomson;
    else
      trace.data = [];
      trace.dim = [];
      trace.dimunits = [];
      trace.x = [];
      trace.t = [];
      trace.units = [];
      trace.name=[num2str(shot) '/' ppftype '/' tracename];
    end
    
  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'te', 'ne'}

    shotfile_exp_eff = AUGexplocation{index};
    
    if strcmp(AUGkeywrdcase{index},'te')
      [a,e]=rdaAUG_eff(shot,'YPR','Te',shotfile_exp_eff);
    else
      [a,e]=rdaAUG_eff(shot,'YPR','Ne',shotfile_exp_eff);
    end
    trace = a;
    trace.data = a.value';
    trace.dim{1} = trace.x;
    trace.dim{2} = trace.t;
    trace.units = trace.unit;
    trace.dimunits{1} = a.area.name ;
    trace.dimunits{2} = a.time_aug.unit;

  case {'equil', 'equil_fpp', 'equil_eqm', 'equil_eqr', 'qrho', 'qrho_fpp'}

    shotfile_exp_eff = AUGexplocation{index};
    
    if strcmp(AUGkeywrdcase{index},'equil_fpp') || strcmp(AUGkeywrdcase{index},'qrho_fpp')
      DIAG = 'FPP';
    elseif strcmp(AUGkeywrdcase{index},'equil_eqm')
      DIAG = 'EQM';
    elseif strcmp(AUGkeywrdcase{index},'equil_eqr')
      DIAG = 'EQR';
    else
      DIAG = 'EQI';
    end
    M_Rmesh_par = sf2par(DIAG,shot,'M','PARMV');
    M_Rmesh = M_Rmesh_par.value + 1; % nb of points
    N_Zmesh_par = sf2par(DIAG,shot,'N','PARMV');
    N_Zmesh = N_Zmesh_par.value + 1; % nb of points
    NTIME_par = sf2par(DIAG,shot,'NTIME','PARMV');
    NTIME = NTIME_par.value; % nb of points
    Lpf_par = rdaAUG_eff(shot,DIAG,'Lpf',shotfile_exp_eff);
    Lpf_tot = Lpf_par.value; % nb of points: 100000*nb_SOL points + nb_core
    Lpf_SOL = fix(Lpf_tot(1:NTIME)/100000);
    Lpf1_t = mod(Lpf_tot(1:NTIME),100000)+1; % nb of Lpf points
% since Lpf depends on time, need to load all first and then loop over time for easier mapping
    [qpsi,e]=rdaAUG_eff(shot,DIAG,'Qpsi',shotfile_exp_eff);
    [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',shotfile_exp);
    [phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',shotfile_exp);
    [Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',shotfile_exp);
    [Area,e]=rdaAUG_eff(shot,DIAG,'Area',shotfile_exp);
    [Ri,e]=rdaAUG_eff(shot,DIAG,'Ri',shotfile_exp);
    [Zj,e]=rdaAUG_eff(shot,DIAG,'Zj',shotfile_exp);
    [PFM_tree,e]=rdaAUG_eff(shot,DIAG,'PFM',shotfile_exp);
    [Pres,e]=rdaAUG_eff(shot,DIAG,'Pres',shotfile_exp);
    [FFP,e]=rdaAUG_eff(shot,DIAG,'FFP',shotfile_exp);
    if isempty(FFP.value); FFP.value=NaN*ones(NTIME,max(Lpf1_t)); end
    [LPFx,e]=rdaAUG_eff(shot,DIAG,'LPFx',shotfile_exp);
    [PFxx,e]=rdaAUG_eff(shot,DIAG,'PFxx',shotfile_exp);
    [RPFx,e]=rdaAUG_eff(shot,DIAG,'RPFx',shotfile_exp);
    [zPFx,e]=rdaAUG_eff(shot,DIAG,'zPFx',shotfile_exp);

    % seems "LCFS" q-value is far too large, limit to some max (when diverted)
    max_qValue = 40; % Note could just put a NaN on LCFS value since ill-defined when diverted
    for it=1:NTIME
      Lpf1 = Lpf1_t(it);
      % Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
      % change it to (radial,time) and use only Lpf+1 points up to LCFS
      if qpsi.value(1,2)<0
	equil.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
      else
	equil.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
      end
      % get x values
      equil.psi(:,it)=psi_tree.value(it,Lpf1:-1:1)';
      equil.psi_axis(it)= equil.psi(1,it);
      equil.psi_lcfs(it)= equil.psi(end,it);
      equil.rhopolnorm(:,it) = sqrt(abs((equil.psi(:,it)-equil.psi_axis(it)) ./(equil.psi_lcfs(it)-equil.psi_axis(it))));
      % get rhotor values
      equil.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)';
      equil.rhotornorm(:,it) = sqrt(abs(equil.phi(:,it) ./ equil.phi(end,it)));
      % get rhovol values
      equil.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)'; 
      equil.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
      equil.rhovolnorm(:,it) = sqrt(abs(equil.vol(:,it) ./ equil.vol(end,it)));
      equil.area(:,it)=Area.value(it,2*Lpf1-1:-2:1)'; 
      equil.dareadpsi(:,it)=Area.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
      equil.Rmesh(:,it) = Ri.value(it,1:M_Rmesh);
      equil.Zmesh(:,it) = Zj.value(it,1:N_Zmesh);
      equil.psi2D(1:M_Rmesh,1:N_Zmesh,it) = PFM_tree.value(1:M_Rmesh,1:N_Zmesh,it);
      equil.pressure(:,it)=Pres.value(it,2*Lpf1-1:-2:1)'; 
      equil.dpressuredpsi(:,it)=Pres.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
      equil.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)';
      equil.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1);
      equil.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1);
      equil.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1);
      %
    end
    equil.x = equil.rhopolnorm;
    %
    [equil_Rcoil,e]=rdaAUG_eff(shot,DIAG,'Rcl',shotfile_exp);
    equil.Rcoils=equil_Rcoil.value;
    [equil_Zcoil,e]=rdaAUG_eff(shot,DIAG,'Zcl',shotfile_exp);
    equil.Zcoils=equil_Zcoil.value;
    % get time values
    [equil_time,e]=rdaAUG_eff(shot,DIAG,'time',shotfile_exp);
    equil.t = equil_time.value(1:NTIME);
    %
    [IpiPSI,e]=rdaAUG_eff(shot,DIAG,'IpiPSI',shotfile_exp);
    equil.Ip = IpiPSI.value';
    %
    equil.data = equil.qvalue; % put q in data
    equil.unit=[]; % not applicable
    equil.name_gdat = [AUGkeywrdcase{index} ' using ' DIAG];
    trace = equil;
    trace.dim{1} = trace.x;
    trace.dim{2} = trace.t;
    trace.units = trace.unit;
    trace.dimunits{1} = '';
    trace.dimunits{2} = 's';

  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'pgyro'}
    %  LOAD MULTI CHANNEL DATA ECS
    %  powers, frequencies, etc 

    shotfile_exp_eff = AUGexplocation{index};

    trace.t=[];
    trace.x=[];
    ppftype='ECS';
    % pgyro tot in index=9
    [a,e]=rdaAUG_eff(shot,ppftype,'PECRH',shotfile_exp_eff);
    if isempty(a) || e~=0
      trace_all = struct([]);
    else
      nb_timepoints = length(a.time_aug.value);
      trace_all.pgyro = NaN*ones(nb_timepoints,9);
      trace_all.freq_ech = NaN*ones(1,8);
      trace_all.pgyro(:,9) = reshape(a.data,nb_timepoints,1);
      trace_all.t = a.time_aug.value;
    end
    for i=1:4
      % "old" ECRH1 gyrotrons: gyro 1 to 4 in pgyro
      tracename_eff = ['PG' num2str(i)];
      [a,e]=rdaAUG_eff(shot,ppftype,tracename_eff,shotfile_exp_eff);
      if isempty(a) || e~=0
      else
	trace_all.ecrh1(i) = a;
	trace_all.pgyro(:,i) = reshape(a.data,nb_timepoints,1);
      end
      a = sf2par('ECS',shot,'gyr_freq',['P_sy1_g' num2str(i)]);
      if isempty(a)
      else
	trace_all.freq_ecrh1(i) = a;
	trace_all.freq_ech(i) = a.value;
      end
      a = sf2par('ECS',shot,'GPolPos',['P_sy1_g' num2str(i)]);
      if isempty(a)
      else
	trace_all.polpos_ecs(i) = a.value;
      end
      a = sf2par('ECS',shot,'GTorPos',['P_sy1_g' num2str(i)]);
      if isempty(a)
      else
	trace_all.torpos_ecs(i) = a.value;
      end
      % "new" ECRH2 gyrotrons: gyro 5 to 8 in pgyro
      tracename_eff = ['PG' num2str(i) 'N'];
      [a,e]=rdaAUG_eff(shot,ppftype,tracename_eff,shotfile_exp_eff);
      if isempty(a) || e~=0
      else
	trace_all.ecrh2(i) = a;
	trace_all.pgyro(:,i+4) = reshape(a.data,nb_timepoints,1);
      end
      a = sf2par('ECS',shot,'gyr_freq',['P_sy2_g' num2str(i)]);
      if isempty(a)
      else
	trace_all.freq_ecrh2(i) = a;
	trace_all.freq_ech(i+4) = a.value;
      end
      a = sf2par('ECS',shot,'GPolPos',['P_sy2_g' num2str(i)]);
      if isempty(a)
      else
	trace_all.polpos_ecs(i+4) = a.value;
      end
      a = sf2par('ECS',shot,'GTorPos',['P_sy2_g' num2str(i)]);
      if isempty(a)
      else
	trace_all.torpos_ecs(i+4) = a.value;
      end
      [a,e]=rdaAUG_eff(shot,'ECN',['G' num2str(i) 'POL'],shotfile_exp_eff);
      if isempty(a) || e~=0
      else
	trace_all.gpol_ecn(i+4) = a;
      end
      [a,e]=rdaAUG_eff(shot,'ECN',['G' num2str(i) 'TOR'],shotfile_exp_eff);
      if isempty(a) || e~=0
      else
	trace_all.gtor_ecn(i+4) = a;
      end
    end
    if ~isempty(trace_all)
      trace_all.dim=[{trace_all.t} {[1:9]}];
      trace_all.data = trace_all.pgyro;
      trace = trace_all;
      trace.x=trace.dim{2};
      trace.dimunits=[{'time [s]'} {'ECRH1(1:4) ECRH2(1:4) ECtot'}];
      trace.units='W';
      trace.freq_ech_units = 'GHz';
      trace.name=[num2str(shot) '/' ppftype '/' 'PGi and PGiN'];
    else
      trace.data = [];
      trace.dim = [];
      trace.dimunits = [];
      trace.x = [];
      trace.t = [];
      trace.units = [];
      trace.freq_ech_units =[]';
     trace.name=[num2str(shot) '/' ppftype '/' 'PGi and PGiN'];
    end
    
  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  case {'powers'}
    %  load powers from TOT timebase as well
    %  

    shotfile_exp_eff = AUGexplocation{index};

    trace.t=[];
    trace.x=[];
    ppftype='TOT';
    % P ohmic
    [a,e]=rdaAUG_eff(shot,ppftype,'P_OH',shotfile_exp_eff);
    if isempty(a) || e~=0
      trace_all = struct([]);
    else
      nb_timepoints = length(a.time_aug.value);
      trace_all.powers = NaN*ones(nb_timepoints,5);
      trace_all.powers(:,1) = reshape(a.data,nb_timepoints,1);
      trace_all.power_names{1} = [ppftype '/P_OH'];
      trace_all.pohmic = a;
      trace_all.t = a.time_aug.value;
    end
    % P NBI
    [a,e]=rdaAUG_eff(shot,ppftype,'PNBI_TOT',shotfile_exp_eff);
    if isempty(a) || e~=0
    else
      trace_all.powers(:,2) = reshape(a.data,nb_timepoints,1);
      trace_all.power_names{2} = [ppftype '/PNBI_TOT'];
      trace_all.pnbi = a;
    end
    % P ECRH
    [a,e]=rdaAUG_eff(shot,ppftype,'PECR_TOT',shotfile_exp_eff);
    if isempty(a) || e~=0
    else
      trace_all.powers(:,3) = reshape(a.data,nb_timepoints,1);
      trace_all.power_names{3} = [ppftype '/PECR_TOT'];
      trace_all.pecrh = a;
    end
    % P ICRH
    [a,e]=rdaAUG_eff(shot,ppftype,'PICR_TOT',shotfile_exp_eff);
    if isempty(a) || e~=0
    else
      trace_all.powers(:,4) = reshape(a.data,nb_timepoints,1);
      trace_all.power_names{4} = [ppftype '/PICR_TOT'];
      trace_all.picrh = a;
    end
    trace_all.powers(:,5) = trace_all.powers(:,1) + trace_all.powers(:,2) + trace_all.powers(:,3) + trace_all.powers(:,4);
    % tau_tot
    [a,e]=rdaAUG_eff(shot,ppftype,'tau_tot',shotfile_exp_eff);
    if isempty(a) || e~=0
    else
      trace_all.tau_tot = a;
      ij=find(~isnan(trace_all.tau_tot.value));
      trace_all.tau_tot_spline=interpos(trace_all.t(ij),trace_all.tau_tot.value(ij),trace_all.t,-1e3);
      trace_all.tau_tot_spline = max(trace_all.tau_tot_spline,0.);
    end
    % betaN
    [a,e]=rdaAUG_eff(shot,ppftype,'beta_N',shotfile_exp_eff);
    if isempty(a) || e~=0
    else
      trace_all.betan = a;
    end
    % betaNthermal
    [a,e]=rdaAUG_eff(shot,ppftype,'beta_Nth',shotfile_exp_eff);
    if isempty(a) || e~=0
    else
      trace_all.betan_thermal = a;
    end
    %
    if ~isempty(trace_all)
      trace_all.dim=[{trace_all.t} {[1:5]}];
      trace_all.data = trace_all.powers;
      trace = trace_all;
      trace.x=trace.dim{2};
      trace.dimunits=[{'time [s]'} {'Pohmic Pnbi Pecrh Picrh Ptot'}];
      trace.units='W';
      trace.name=[num2str(shot) '/' ppftype '/' 'P_OH, PNBI_TOT, PECR_TOT, PICR_TOT'];
    else
      trace.data = [];
      trace.dim = [];
      trace.dimunits = [];
      trace.x = [];
      trace.t = [];
      trace.units = [];
      trace.freq_ech_units =[]';
     trace.name=[num2str(shot) '/' ppftype '/' 'P_OH, PNBI_TOT, PECR_TOT, PICR_TOT'];
    end
    
  otherwise
    disp('case not yet defined')
    
end
