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}!)
% 'SXB'  =  soft x-ray emission from (by default camera J) SXB/J_xx camera (sxb, sxB, etc all work)
%
%     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', 'SXB/J_053', 'SXB'
%           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)
%                  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');
%

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,[{'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(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
  if ~isempty(strmatch(data_type_eff_noext,[{'Ha'} {'ha'} {'Halpha'} {'halpha'}],'exact'))
    data_type_eff_noext='Halpha';
  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'} {'qrho_cliste'} {'q95'} {'kappa'} ...
      {'delta'} {'deltatop'} {'deltabot'} {'neint'} ...
      {'ne'} {'te'} {'nerho'} {'terho'} ...
		     {'sxr'} {'sxR'} {'sxb'} {'ece'} {'Halpha'}];
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.iqrho_cliste=strmatch('qrho_cliste',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.isxb=strmatch('sxb',AUGkeywrdall,'exact');
AUGsig.iece=strmatch('ece',AUGkeywrdall,'exact');
AUGsig.iHalpha=strmatch('Halpha',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_cliste)=AUGkeywrdall(AUGsig.iqrho_cliste); % 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.isxb)=AUGkeywrdall(AUGsig.isxb);
%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'};
AUGsiglocation(:,AUGsig.iece)={'CEC'; 'Trad-A'};
AUGsiglocation(:,AUGsig.iHalpha)={'POT'; 'ELMa-Han'};

% 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
    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,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 
    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,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 {'te', 'ne'}

    if strcmp(AUGkeywrdcase{index},'te')
      [a,e]=rdaAUG_eff(shot,'YPR','Te');
    else
      [a,e]=rdaAUG_eff(shot,'YPR','Ne');
    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 {'qrho', 'qrho_cliste'}

    if strcmp(AUGkeywrdcase{index},'qrho')
      DIAG = 'FPP';
    else
      DIAG = 'EQI';
    end
    [a,e]=rdaAUG_eff(shot,DIAG,'Qpsi');
    % Qpsi has inverted channel/time from CEC
    a.value = a.value(:,end:-1:1)';
    a.data = a.value;
    % get x values
    [psi,e]=rdaAUG_eff(shot,DIAG,'PFL');
    psi.value=psi.value(:,end:-1:1)';
    psi_axis= (ones(size(psi.value,1),1) * psi.value(1,:));
    a.x = sqrt(abs((psi.value-psi_axis) ./(zeros(size(psi.value))-psi_axis) )); % psi_edge=0 assumed
    a.psi = psi.value;
    % get time values
    [psi_time,e]=rdaAUG_eff(shot,DIAG,'time');
    if length(psi_time.value) > length(a.x(1,:))
      a.t = psi_time.value(1:length(a.x(1,:))); % problem with times??
    elseif length(psi_time.value) < length(a.x(1,:))
      a.t = psi_time.value; % problem with times??
      a.t(end+1:length(a.x(1,:))) = psi_time.value(end) + [1:length(a.x(1,:))-length(psi_time.value)].*diff(psi_time.value(end-1:end));
    else
      a.t = psi_time.value;
    end
    % get rhotor values
    [phi,e]=rdaAUG_eff(shot,DIAG,'TFLx');
    phi.value=phi.value(:,end:-1:1)';
    a.rhotor = sqrt(abs(phi.value ./ (ones(size(phi.value,1),1) * phi.value(end,:))));
    a.torflux = phi.value;
    % get rhovol values
    [Vol,e]=rdaAUG_eff(shot,DIAG,'Vol');
    Vol.value=Vol.value(:,end-1:-2:1)'; % 2nd index are dV/dpsi
    a.rhovol = sqrt(abs(Vol.value ./ (ones(size(Vol.value,1),1) * Vol.value(end,:))));
    a.volume = Vol.value;
    %
    trace = a;
    trace.dim{1} = trace.x;
    trace.dim{2} = trace.t;
    trace.units = trace.unit;
    trace.dimunits{1} = '';
    trace.dimunits{2} = 's';

  otherwise
    disp('case not yet defined')
    
end