rdaAUG_eff.m

function [adata,error]=rdaAUG_eff(shot,diagname,sigtype,shotfile_exp,varargin);
%
% [adata,error]=rdaAUG_eff(shot,diagname,sigtype,shotfile_exp,varargin);
%
% inputs: shot,diagname,sigtype,shotfile_exp are mandatory, so specify default shotfile_exp before for example
%
% gets data using sf2sig or mdsplus (when mdsplus will be available)
% 1D arrays: assumes dimension is time
% 2D arrays: assumes data vs (x,time)
% 3D arrays: assumes data vs (x,time,hsig) (for mdsplus)
%
% varargin{1}: time interval or timevalue, will get data closest to that time or within that time interval
% varargin{2}: extra args for sf2sig like '-raw'
%              if starts as in varargin{3}, assume varargin{2} is actually varargin{3} (allows param:... set in gdat easily)
%
% varargin{3}: 'param:xxx' means get parameter (instead of signal) with xxx as parameter to get
%              'param-set:xxx' means get parameter-set (at this stage available only with sf2ps)
%              'area-base:yy:i' means it's an area-base type signal and if known from yy dim_of(i)
%              'time-base:yy:i' means it's a  time-base type signal and if known from yy dim_of(i)
%
%     for mds usage, calls: augparam (_shot, _diag, _psetname, _parameter, _experiment, _edition, _oshot, _oedition)
%     for sf2 usage: help sf2par or sf2ps on ipp relevant node
%            : [] empty means get signal (default),
%     for mds usage: augdiag (_shot, _diag, _signame, _experiment, _edition, _t1, _t2, _oshot, _oedition)
%     for sf2 usage: help sf2sig on ipp relevant node
%
% examples:
%          [data,error]=rdaAUG_eff(15133,'MAG','Ipi');
%          [data,error]=rdaAUG_eff(15133,'MAG','Ipi',[1 5]);
%
% set global variable: usemdsplus to decide if sf2sig or mdsplus is used:
%     usemdsplus=1 except if sf2sig is a function, then usemdsplus=0
%    (----begin of old comments: >> global usemdsplus
%        >> usemdsplus=1 % means use mds to get data
%        >> usemdsplus=0 % means use sf2sig  (default if not defined)
%       if ~exist('usemdsplus'); usemdsplus=0; end   -----end of old comments)
%
%
% return: 1D array and dimensions as 1xN
%         2D array transposed (assumed were changed from C origin, this way gets usually time as 2nd)
%
%global usemdsplus
usemdsplus = ~[exist('sf2sig')==3];
if ~exist('usemdsplus') || isempty(usemdsplus); usemdsplus=0; end
error=999;
adata.data = [];
adata.value = [];
adata.units = [];
adata.dim = [];
adata.dimunits = [];
adata.t = [];
adata.x = [];
adata_time.data = [];
adata_time.value = [];
adata_area = [];
adata.time_aug = adata_time;
adata.area = adata_area;
adata.exp = shotfile_exp;
adata.edition_in = [];
adata.edition_out = [];

varargin_eff = varargin;
nargin_eff = nargin;
% now use special_signal (4th input) to specify extra demand, so can remain with edition or raw in extra_arg
% $$$ special_varargin3 = {'param:', 'param-set:', 'area-base:', 'time-base:'};
% $$$ ab=regexpi(varargin_eff{2},special_varargin3);
% $$$ if ~isempty(cell2mat(ab)) && (length(varargin_eff)<3 || isempty(varargin_eff{3}))
% $$$   varargin_eff{3} = varargin_eff{2};
% $$$   varargin_eff{2} = '';
% $$$   nargin_eff = max(nargin,7);
% $$$ end

time_int=[];
if nargin_eff>=5 & ~isempty(varargin_eff{1})
  time_int=varargin_eff{1};
end
extra_arg_sf2sig=[];
if nargin_eff>=6 & ~isempty(varargin_eff{2})
  extra_arg_sf2sig=varargin_eff{2};
end
param_name=[];
param_set_name=[];
area_base=false;
area_base_name=[];
area_base_dimof=[];
time_base=false;
time_base_name=[];
time_base_dimof=[];
if nargin_eff>=7 && ~isempty(varargin_eff{3}) && ischar(varargin_eff{3}) ...
           && length(varargin_eff{3})>=7 && strcmp(lower(varargin_eff{3}(1:6)),'param:')
  param_name=varargin_eff{3}(7:end);
end
if nargin_eff>=7 && ~isempty(varargin_eff{3}) && ischar(varargin_eff{3}) ...
           && length(varargin_eff{3})>=11 && strcmp(lower(varargin_eff{3}(1:10)),'param-set:')
  param_set_name=varargin_eff{3}(11:end);
end
if nargin_eff>=7 && ~isempty(varargin_eff{3}) && ischar(varargin_eff{3}) ...
           && length(varargin_eff{3})>=9 && strcmp(lower(varargin_eff{3}(1:9)),'area-base')
  area_base = true;
  ij=findstr(varargin_eff{3},':');
  if length(ij)==2
    area_base_name=varargin_eff{3}(ij(1)+1:ij(2)-1);
    area_base_dimof = str2num(varargin_eff{3}(ij(2)+1:end));
  elseif length(ij)==1
    area_base_name=varargin_eff{3}(ij(1)+1:end);
  end
end
if nargin_eff>=7 && ~isempty(varargin_eff{3}) && ischar(varargin_eff{3}) ...
           && length(varargin_eff{3})>=9 && strcmp(lower(varargin_eff{3}(1:9)),'time-base')
  time_base = true;
  ij=findstr(varargin_eff{3},':');
  if length(ij)==2
    time_base_name=varargin_eff{3}(ij(1)+1:ij(2)-1);
    time_base_dimof = str2num(varargin_eff{3}(ij(2)+1:end));
  elseif length(ij)==1
    time_base_name=varargin_eff{3}(ij(1)+1:end);
  end
end

did_transpose = 0;
if usemdsplus
  % a_remote=mdsremotelist; does not seem sufficient if connection broken, test with 1+2
  is3 = mdsvalue('1+2');
  if is3 ~= 3
    warning('not connected to an mds server, cannot use mds to get data')
    return
  end
  % use mdsplus

% $$$   if ~unix('test -d /home/duval/mdsplus')
% $$$     addpath('/home/duval/mdsplus')
% $$$   end

% $$$   if ~unix('test -d /home/osauter/gdat')
% $$$     mdsconnect('localhost:8001');
% $$$   else
% $$$     mdsconnect('localhost');
% $$$   end

  % extract if need raw
  ij=strfind(extra_arg_sf2sig,'raw');
  ask_raw = '';
  if ~isempty(ij)
    ask_raw = '"raw"';
  end
  % extract edition number if provided as '-ed',value in extra_arg_sf2sig
  ij=strfind(extra_arg_sf2sig,'-ed');
  ed_number = '';
  if ~isempty(ij)
    ed_number = num2str(sscanf(extra_arg_sf2sig(ij+5:end),'%d'));
  end
  if isempty(time_int)
    tstart = '';
    tend = '';
  else
    tstart = num2str(time_int(1));
    tend = num2str(time_int(2));
  end
  user=getenv('USER');
  if isempty(param_name) && isempty(param_set_name) && ~area_base && ~time_base
    % use augsignal to get effective layout as in ISIS and sf2sig, for example for EQI/PFM
    eval(['[data,error]=mdsvalue(''_rdaeff' user diagname '=augsignal(' num2str(shot) ',"' diagname '","' sigtype '","' shotfile_exp ...
            '",' ed_number ',' tstart ',' tend ',_oshot' user diagname ',_oed' user diagname ',' ask_raw ')'');']);
  elseif isempty(param_set_name) && ~area_base && ~time_base
    % use augparam
    eval(['[data,error]=mdsvalue(''_rdaeff' user diagname '=augparam(' num2str(shot) ',"' diagname '","' sigtype '","' param_name '"' ...
          ',"' shotfile_exp '",' ed_number ',_oshot' user diagname ',_oed' user diagname ')'');']);
  elseif ~area_base && ~time_base
    % param-set, cannot get this yet with mdsvalue
    disp(['cannot get param-set with mds yet (only sf2ps): ' param_set_name])
    data = [];
    error = 11;
  else
    % area-base or time-base, can only get dim_of at this stage
    area_time_base_name = area_base_name;
    area_time_base_dimof = area_base_dimof;
    if ~isempty(time_base_name);
      area_time_base_name = time_base_name;
      area_time_base_dimof = time_base_dimof;
    end
    if ~isempty(area_time_base_name)
      if nargin_eff>=5 & ~isempty(varargin_eff{1})
        eval(['[data,error]=mdsvalue(''_rdaeff' user diagname '=augdiag(' num2str(shot) ',"' diagname '","' area_time_base_name '","' shotfile_exp ...
              '",' ed_number ',' num2str(varargin_eff{1}(1),'%.14f') ',' num2str(varargin_eff{1}(end),'%.14f') ')'');']);
      else
        eval(['[data,error]=mdsvalue(''_rdaeff' user diagname '=augdiag(' num2str(shot) ',"' diagname '","' area_time_base_name '","' shotfile_exp ...
              '",' ed_number ',,,_oshot' user diagname ',_oed' user diagname ')'');']);
      end
      if ~isempty(area_time_base_dimof)
        eval(['data=mdsvalue(''dim_of(_rdaeff' user diagname ',' num2str(area_time_base_dimof) ')'');']);
      else
        for j=1:length(size(data))
          eval(['dataj=mdsvalue(''dim_of(_rdaeff' user diagname ',' num2str(j) ')'');']);
          if (prod(size(dataj))~=length(dataj))
            data = dataj;
            area_time_base_dimof = j;
            break
          end
        end
      end
    end

  end
  adata.data=data;
  if exist('ed_number')
    adata.edition_in = str2num(ed_number);
  end
  adata.edition_out = mdsvalue(['_oed' user diagname]);
  if mod(error,2)==2; disp(['error even for ' diagname ' ; ' sigtype]); end
  % begin some fix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % fix some problems (seems not needed anymore with lxmdsplus.aug.ipp.mpg.de, but do with mdsplus.aug.ipp.mpg.de (so leave for some timing tests)
  nb_surfaces_m1 = 40;
  nb_xpoints_m1 = 4;
  if strcmp(sigtype,'Lpf')
    if (min(adata.data)<0*100000+3 && min(adata.data)~=0) || max(adata.data)>100000*20
      disp(['seems to be a problem with Lpf assume 1 sol point and ' num2str(nb_surfaces_m1+1) ' surfaces'])
      adata.data(:) = 100000 + nb_surfaces_m1;
    end
  end
  if strcmp(sigtype,'LPFx')
    if (min(adata.data)<1 && min(adata.data)~=0) || max(adata.data)>1000
      disp(['seems to be a problem with LPFx assume ' num2str(nb_xpoints_m1+1) ' special X-points'])
      adata.data(:) = nb_xpoints_m1;
    end
  end
  % end of some fix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  hsig=[];
  ss=size(data);
  nbofdim=length(ss);
  if prod(ss)==length(data); nbofdim=1; end
  nbofdim=max(nbofdim,1);
  switch nbofdim
   case 1
    adata.data=reshape(adata.data,1,length(adata.data));
    idim0 = 0;
    eval(['time=mdsvalue(''dim_of(_rdaeff' user diagname ',' num2str(idim0) ')'');']);
    if numel(time) ~= numel(adata.data)
      idim0 = 1;
      eval(['time=mdsvalue(''dim_of(_rdaeff' user diagname ',' num2str(idim0) ')'');']);
      if numel(time) ~= numel(adata.data) && (ischar(adata.data) && isempty(strfind(lower(adata.data),'abort')))
        warning(['problem with dim for: ' diagname ', ' sigtype])
      end
    end
    time = reshape(time,1,length(time));
    x=[];
    adata.dim = {time};
    eval(['tunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',' num2str(idim0) '))'');']);
    adata.dimunits = {tunits};
   case 2
    if strcmp(upper(diagname),'IDA')
      % rho, time in original dimension, so do not transpose
      idim_x = 0;
      idim_t = 1;
    else
      adata.data = adata.data';
      did_transpose = 1;
      idim_x = 1;
      idim_t = 0;
    end
    % transposed because of C relation and backward compatibility with sf2sig part
    eval(['x=mdsvalue(''dim_of(_rdaeff' user diagname ',' num2str(idim_x) ')'');']);
    if prod(size(x))==length(x); x = reshape(x,1,length(x)); end
    eval(['time=mdsvalue(''dim_of(_rdaeff' user diagname ',' num2str(idim_t) ')'');']);
    time = reshape(time,1,length(time));
    adata.dim = {x, time};
    eval(['xunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',' num2str(idim_x) '))'');']);
    eval(['tunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',' num2str(idim_t) '))'');']);
    adata.dimunits = {xunits, tunits};

% $$$    case 3
% $$$     eval(['x=mdsvalue(''dim_of(_rdaeff' user diagname ',0)'');']);
% $$$     if prod(size(x))==length(x); x = reshape(x,1,length(x)); end
% $$$     eval(['time=mdsvalue(''dim_of(_rdaeff' user diagname ',1)'');']);
% $$$     time = reshape(time,1,length(time));
% $$$     disp('3rd dimension in hsig!!!!!!!!!!!!!!!!!!!!!!!!!')
% $$$     eval(['hsig=mdsvalue(''dim_of(_rdaeff' user diagname ',2)'');']);
% $$$     if prod(size(hsig))==length(hsig); hsig = reshape(hsig,1,length(hsig)); end
% $$$     adata.dim = {x, time, hsig};
% $$$     eval(['xunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',0))'');']);
% $$$     eval(['tunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',1))'');']);
% $$$     eval(['hsigunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',2))'');']);
% $$$     adata.dimunits = {xunits, tunits, hsigunits};
% $$$     [zz,itime] = max(size(adata.data));
% $$$     for i=1:nbofdim
% $$$       if strcmp(adata.dimunits{i},'s'); itime = i; end
% $$$     end
% $$$     ix = min(setdiff([1:2],itime));
% $$$     ihsig = setdiff([1:nbofdim],[ix itime])
% $$$     x = adata.dim{ix};
% $$$     time = adata.dim{itime};
% $$$     hsig = adata.dim{ihsig};
% $$$     adata.dim = {x, time, hsig};
% $$$     xunits = adata.dimunits{ix};
% $$$     tunits = adata.dimunits{itime};
% $$$     hsigunits = adata.dimunits{ihsig};
% $$$     adata.dimunits = {xunits, tunits, hsigunits};

   otherwise
    itime = 1; % default
    for i=1:nbofdim
      eval(['dimarray=mdsvalue(''dim_of(_rdaeff' user diagname ',' num2str(i-1) ')'');']);
      if prod(size(dimarray)) == length(dimarray)
        eval(['adata.dim{' num2str(i) '}=reshape(dimarray,1,length(dimarray));']);
      else
          eval(['adata.dim{' num2str(i) '} = dimarray;']);
      end
      eval(['adata.dimunits{' num2str(i) '}=deblank(mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',' num2str(i-1) '))''));']);
      if strcmp(adata.dimunits{i},'s'); itime = i; end
    end
    x = adata.dim{min(setdiff([1:2],itime))};
    time = adata.dim{itime};
  end

  adata.value = adata.data; % for backward compatibility might need to keep .value, to be checked later, .data for sure
  adata.t=time;
  adata.x=x;
  adata.hsig=hsig;
  eval(['adata.units=deblank(mdsvalue(''units_of(_rdaeff' user diagname ')''));']);
  %mdsdisconnect;
% $$$   if ~unix('test -d /home/duval/mdsplus')
% $$$     rmpath('/home/duval/mdsplus')
% $$$   end

else
  % use sf2sig, sf2par, sf2ab or sf2tb
  adata_time.data = [];
  adata_area = [];
  if isempty(param_name) && isempty(param_set_name) && ~area_base && ~time_base
    % use sf2sig
    if isempty(time_int)
      try
        if ~isempty(extra_arg_sf2sig) && ~strcmp(extra_arg_sf2sig,'[]')
          eval(['[adata,adata_time, adata_area]=sf2sig(diagname,shot,sigtype,''-exp'',shotfile_exp,' extra_arg_sf2sig ');']);
        else
          [adata,adata_time, adata_area]=sf2sig(diagname,shot,sigtype,'-exp',shotfile_exp);
        end
	adata.data = adata.value; % at this stage keep both but try to work only on .data
      catch ME
        throw(ME)
      end
    else
      try
        if ~isempty(extra_arg_sf2sig) && ~strcmp(extra_arg_sf2sig,'[]')
          eval(['[adata,adata_time, adata_area]=sf2sig(diagname,shot,sigtype,[' num2str(time_int(1)) ';' num2str(time_int(2)) ...
                '] ,''-exp'',shotfile_exp,' extra_arg_sf2sig ');']);
        else
          [adata,adata_time, adata_area]=sf2sig(diagname,shot,sigtype,[time_int(1);time_int(end)],'-exp',shotfile_exp);
        end
          adata.data = adata.value; % at this stage keep both but try to work only on .data (for adata which can be from sf2.. or mds)
      catch ME
        throw(ME)
      end
    end

    if isempty(adata.data)
      return
    end

    %%%%%%%%%%%%%%%%%%%%%%%% NEED TO DO THESE AFTER EITHER MDSPLUS or SF2SIG, after data and dims obtained %%%%%%%%%%%%%%%%%%%%%%%
    % special checks
    if strcmp(upper(diagname),'SXB')
      % time missing one point
      if length(adata.data) == length(adata_time.value)+1
        adata_time.value=linspace(adata_time.range(1),adata_time.range(2),length(adata.data));
        adata_time.index(2) = length(adata.data);
      end
    end
    %
% $$$     if strcmp(upper(sigtype),'PNIQ')
% $$$       % transform 4x2 PINIs in 1:8 PINIs and total in index=9
% $$$       if (prod(size(adata.data))/length(adata_time.value) == 8)
% $$$         tmp(:,1:4) = adata.data(:,:,1);
% $$$         tmp(:,5:8) = adata.data(:,:,2);
% $$$         tmp(:,9) = sum(tmp,2);
% $$$         adata.data = tmp'; % transpose since will be transposed afterwards
% $$$         adata.dimunits = {'s','8 sources;total'};
% $$$       else
% $$$         disp('expects 8 sources in PNIQ');
% $$$         return
% $$$       end
% $$$     end

    adata.time_aug = adata_time;

    adata.area = adata_area;

    adata.exp = shotfile_exp;
    if (prod(size(adata.data))==length(adata.data))
      % only time signal
      adata.x = [];
      adata.data=reshape(adata.data,1,length(adata.data));
      if ~isempty(adata.time_aug)
        adata.t=adata.time_aug.value;
      else
        adata.t=[1:size(adata.data,2)];
      end
    else
      if strcmp(upper(diagname),'IDA')
        % rho, time in original dimension, so do not transpose
      else
        if length(size(adata.data))<=2; 
          adata.data = adata.data'; 
          did_transpose = 1;
        end % cannot transpose Nd>2 matrix
      end
      if ~isempty(adata.time_aug)
        if length(size(adata.data))<=2;
          adata.x=[1:prod(size(adata.data))/length(adata_time.value)];
        else
          adata.x = [];
        end
        adata.t=adata.time_aug.value;
      else
        adata.x=[1:size(adata.data,1)];
        adata.t=[1:size(adata.data,2)];
      end
      if ~isempty(adata.area)
        adata.x = adata.area.value{1};
      end
    end
    adata.units = adata.unit;
  elseif isempty(param_set_name) && ~area_base && ~time_base
    % use sf2par
    try
      if ~isempty(extra_arg_sf2sig) && ~strcmp(extra_arg_sf2sig,'[]')
	eval(['[adata]=sf2par(diagname,shot,param_name,sigtype,''-exp'',shotfile_exp,' extra_arg_sf2sig ');']);
      else
	eval(['[adata]=sf2par(diagname,shot,param_name,sigtype,''-exp'',shotfile_exp);']);
      end
      adata.data = adata.value;
    catch ME
      throw(ME)
    end
  elseif ~area_base && ~time_base
    % use sf2ps
    try
      if ~isempty(extra_arg_sf2sig) && ~strcmp(extra_arg_sf2sig,'[]')
	eval(['[adata]=sf2ps(diagname,shot,param_set_name,sigtype,''-exp'',shotfile_exp,' extra_arg_sf2sig ');']);
      else
	eval(['[adata]=sf2ps(diagname,shot,param_set_name,sigtype,''-exp'',shotfile_exp);']);
      end
      adata.data = adata.value;
    catch ME
      throw(ME)
    end
  elseif ~time_base
    % area-base
    try
      if ~isempty(extra_arg_sf2sig) && ~strcmp(extra_arg_sf2sig,'[]')
        eval(['[adata]=sf2ab(diagname,shot,area_base_name,''-exp'',shotfile_exp,' extra_arg_sf2sig ');']);
      else
        [adata]=sf2ab(diagname,shot,area_base_name,'-exp',shotfile_exp);
      end
      adata.data = adata.value{1};
      adata.value = adata.data;
    catch ME
      throw(ME)
    end
  else
    % time-base
    try
      if ~isempty(extra_arg_sf2sig) && ~strcmp(extra_arg_sf2sig,'[]')
        eval(['[adata]=sf2tb(diagname,shot,time_base_name,''-exp'',shotfile_exp,' extra_arg_sf2sig ');']);
      else
        [adata]=sf2tb(diagname,shot,time_base_name,'-exp',shotfile_exp);
      end
      adata.data = adata.value;
    catch ME
      throw(ME)
    end
  end
end

adata.value = adata.data;

if strcmp(upper(sigtype),'PNIQ')
  % transform 4x2 PINIs in 1:8 PINIs and total in index=9
  if (prod(size(adata.data))/length(adata.t) == 8)
    tmp(:,1:4) = adata.data(:,:,1);
    tmp(:,5:8) = adata.data(:,:,2);
    tmp(:,9) = sum(tmp,2);
    adata.data = tmp;
    adata.value = adata.data;
    adata.x = [1:9];
    adata.dim = {adata.t, adata.x};
    adata.dimunits = {'s','8 sources;total'};
  else
    disp('expects 8 sources in PNIQ');
    return
  end
end

if strcmp(diagname,'IDA') && strcmp(sigtype,'rhot')
  % In some shots, rhot was multiplied by 1e6, so fix here:
  if max(max(adata.data)) > 1e5
    adata.data = adata.data / 1e6;
  end
end

error=0;