loadAUGdata.m

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}!)
%
%     gdat(15133,'MAG/Ipi',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'
%           data_type='POT/ELMa-Han', 'MOD/OddNAmp', 'MOD/EvenNAmp', 'TOT/PNBI_TOT'
%
% Meaning 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
%                  varargin{4}:  time range [t1 t2] (to limit data collected)
%
% 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');
%

varargout={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;
if size(data_type_eff,1)==1
  i=findstr('/',data_type_eff);
  if length(i)>=1
    % assumes given a la 'MAG/Ipi'
    data_type_eff=[{data_type_eff(1:i(1)-1)} ; {data_type_eff(i(1)+1:end)}];
  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'} {'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,[{'Terho'}],'exact'))
    data_type_eff_noext='terho';
  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,[{'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,[{'Rcont'}],'exact'))
    data_type_eff_noext='rcont';
  end
  if ~isempty(strmatch(data_type_eff_noext,[{'Zcont'}],'exact'))
    data_type_eff_noext='zcont';
  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'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'qrho'} {'q95'} {'kappa'} ...
      {'delta'} {'deltatop'} {'deltabot'} {'neint'} ...
      {'ne'} {'te'} {'nerho'} {'terho'} ...
      {'sxr'} {'sxR'} {'ece'}];
AUGsig.iip=strmatch('Ip',AUGkeywrdall,'exact');
AUGsig.izmag=strmatch('zmag',AUGkeywrdall,'exact');
AUGsig.irmag=strmatch('rmag',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.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.isxr=strmatch('sxr',AUGkeywrdall,'exact');
AUGsig.isxR=strmatch('sxR',AUGkeywrdall,'exact');
AUGsig.iece=strmatch('ece',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.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.iterho)=AUGkeywrdall(AUGsig.iterho); % idem
AUGkeywrdcase(AUGsig.isxr)=AUGkeywrdall(AUGsig.isxr);
AUGkeywrdcase(AUGsig.isxR)=AUGkeywrdall(AUGsig.isxR);
%AUGkeywrdcase(AUGsig.iece)=AUGkeywrdall(AUGsig.iece);

% 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'; 'Ipi'};
AUGsiglocation(:,AUGsig.izmag)={'FPG'; 'Zmag'};
AUGsiglocation(:,AUGsig.irmag)={'FPG'; 'Rmag'};
AUGsiglocation(:,AUGsig.ircont)={'' ; ''}; AUGsigtimeindx(AUGsig.ircont)=2;
AUGsiglocation(:,AUGsig.izcont)={'' ; ''}; AUGsigtimeindx(AUGsig.izcont)=2;
AUGsiglocation(:,AUGsig.ivol)={''; ''};
AUGsiglocation(:,AUGsig.iq95)={''; ''};
AUGsiglocation(:,AUGsig.ikappa)={''; ''};
AUGsiglocation(:,AUGsig.ideltatop)={''; ''};
AUGsiglocation(:,AUGsig.ideltabot)={''; ''};
AUGsiglocation(:,AUGsig.ineint)={'DCN'; 'H-1'};

% 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}];
    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};
    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);
%    switch tracename
%  special cases if traces do not exist for some shot or other
%    end
    
    trace=a;
    clear error
    error=e;
    if 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 
    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
    trace.t=[];
    trace.x=[];
    ppftype='SXR';
    tracename='B';
    [a,e]=rdaAUG_eff(shot,ppftype,tracename,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
    nth=13;
    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);
      [R_B, Z_B, ang_B,Rsxr]=sxrbgeometry(zmageff);
      radius.data=Rsxr;
      radius.t=trace.t;
      varargout{1}={radius};
      trace.R=radius;
    end

  otherwise
    disp('case not yet defined')
    
end