Check cxrs

I removed CXRS files to use cxrs git now, needs to be merged with new master now

I'll reopen a branch if needed to get other stuff from cxrs for conf, for Topic 01 and IDS

matlab/TCV/gdat_tcv.m

@@ -1939,11 +1939,11 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ne','te'}
     % ne or Te from Thomson data on raw z mesh vs (z,t)
-    if ~(isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
-         gdat_data.gdat_params.edge>0)
-      gdat_data.gdat_params.edge = 0;
+    if ~(isfield(gdat_data.gdat_params,'systems') && ~isempty(gdat_data.gdat_params.systems) && ...
+        all(ismember(gdat_data.gdat_params.systems,{'all','main','edge'})))
+      gdat_data.gdat_params.systems = {'all'};
     end
-    [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.edge,gdat_params.nverbose);
+    [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.systems,gdat_params.nverbose);
     if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
       [gdat_data] = gdat2time_out(gdat_data,1,{'error_bar'});
     end
@@ -1957,35 +1957,31 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     else
       gdat_data.gdat_params.zshift = zshift;
     end
-    if ~(isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
-         gdat_data.gdat_params.edge>0)
-      gdat_data.gdat_params.edge = 0;
+    if ~(isfield(gdat_data.gdat_params,'systems') && ~isempty(gdat_data.gdat_params.systems) && ...
+        all(ismember(gdat_data.gdat_params.systems,{'all','main','edge'})))
+      gdat_data.gdat_params.systems = {'all'};
     end
-    [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.edge,gdat_params.nverbose);
+    [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.systems,gdat_params.nverbose);
     % construct rho mesh
-    edge_str_ = '';
-    edge_str_dot = '';
-    if gdat_data.gdat_params.edge
-      edge_str_ = '_edge';
-      edge_str_dot = '.edge';
-    end
-    if liuqe_matlab==1 && strcmp(substr_liuqe,'_1'); substr_liuqe = ''; end
     % psiscatvol obtained from linear interpolation in fact so not quite ok near axis where psi is quadratic
     % NOTE: Since 2019 psi_scatvol uses psitbx with cubic-spline interpolation, unsure what shots have new data or not ...
     recompute_psiscatvol_always = 1;
     if liuqe_version==-1; recompute_psiscatvol_always = 1; end
     % NOTE: Since oct.2019 psi_scatvol stored in nodes now uses LIUQE.M (shot range to check if any ...)
     if all(abs(zshift)<1e-5) && liuqe_matlab==0 && recompute_psiscatvol_always==0
-      psi_max=gdat_tcv([],['\results::thomson' edge_str_dot ':psi_max' substr_liuqe],'nverbose',gdat_params.nverbose);
-      psiscatvol=gdat_tcv([],['\results::thomson' edge_str_dot ':psiscatvol' substr_liuqe],'nverbose',gdat_params.nverbose);
+      if strcmp(substr_liuqe,'_1'); substr_liuqe = ''; end
+      [psiscatvol,psi_max] = get_thomson_psiscatvol(shot,gdat_data.gdat_params.systems,substr_liuqe,gdat_params.nverbose);
+      % NOTE: time and Z are switched here (same as previous state) and
+      % code is not reachable anyway) ... waiting for a fix
     else
       % calculate new psiscatvol
-      [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, edge_str_dot, psitbx_str, gdat_params);
-
-      % Add the results to the output of gdat
-      gdat_data.psiscatvol = psiscatvol;
-      gdat_data.psi_max = psi_max;
+      [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, psitbx_str, gdat_params);
     end
+    
+    % Add the results to the output of gdat
+    gdat_data.psiscatvol = psiscatvol;
+    gdat_data.psi_max = psi_max;
+    
     if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data)
       psi_norm = 1.-psiscatvol.data./repmat(psi_max.data(:).',[size(psiscatvol.data,1),1]);
       mask = psi_norm < 0;
@@ -2020,15 +2016,12 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       gdat_data.ne.units = 'm^{-3}';
       gdat_data = rmfield(gdat_data,{'firrat','data_raw','error_bar_raw'});
       %
-      nodenameeff=['\results::thomson' edge_str_dot ':te'];
-      tracetdi=tdi(nodenameeff);
-      nodenameeff=['\results::thomson' edge_str_dot ':te; error_bar'];
-      tracestd=tdi(['\results::thomson' edge_str_dot ':te:error_bar']);
-      gdat_data.te.data=tracetdi.data';
-      gdat_data.te.error_bar=tracestd.data';
-      gdat_data.te.units = tracetdi.units;
-      gdat_data.data_fullpath=['pe=1.6e-19*ne*Te in data, .ne, .te from \results::thomson' ...
-                    edge_str_dot ':ne and te and projected on rhopol\_norm'];
+      te = get_thomson_raw_data(shot,'te',gdat_data,gdat_data.gdat_params.systems,gdat_params.nverbose);
+      gdat_data.te.data=te.data;
+      gdat_data.te.error_bar=te.error_bar;
+      gdat_data.te.units = te.units;
+      gdat_data.data_fullpath=['pe=1.6e-19*ne*Te in data, .ne, .te from \results::thomson for systems (' ...
+                    strjoin(gdat_data.gdat_params.systems,'/') ') and projected on rhopol\_norm'];
       gdat_data.units='N/m^2; 1.6022e-19 ne Te';
       gdat_data.data = 1.6022e-19 .* gdat_data.ne.data .* gdat_data.te.data;
       gdat_data.error_bar = 1.6022e-19 .* (gdat_data.ne.data .* gdat_data.te.error_bar ...
@@ -2077,7 +2070,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         def_rhotornorm = '\results::conf:rhotor';
         def_rhovolnorm = '\results::conf:rhovol';
        otherwise
-        if (gdat_params.nverbose>=1);
+        if (gdat_params.nverbose>=1)
           disp('should not be in switch gdat_data.gdat_params.fit_type')
           disp('ask olivier.sauter@epfl.ch')
         end
@@ -2091,7 +2084,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         elseif strcmp(data_request_eff(1:2),'te')
           nodenameeff = [def_proffit 'teft'];
         else
-          if (gdat_params.nverbose>=1);
+          if (gdat_params.nverbose>=1)
             disp(['should not be here: data_request_eff, data_request_eff(1:2)= ',data_request_eff, data_request_eff(1:2)]);
           end
         end
@@ -2506,7 +2499,6 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       % add extra fields then correct
       [gdat_data] = gdat2time_out(gdat_data,21);
     end
-
    case {'phi_tor', 'phitor', 'toroidal_flux'}
     % Phi(LCFS) = int(Bphi dSphi), can use Eq.(11) of "Tokamak coordinate conventions: COCOS" paper:
     % O. Sauter, S.Yu. Medvedev, Comput. Phys. Commun. 184 (2013) 293–302
@@ -2635,7 +2627,6 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
       [gdat_data] = gdat2time_out(gdat_data,1);
     end
-
    case {'pressure', 'pressure_rho', 'p_rho'}
     if liuqe_matlab==0
       nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('p_rho',1,0) '",''' psitbx_str ''')'];
@@ -3397,62 +3388,38 @@ error_status=0;
 return
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function [firthomratio] = get_fir_thom_rat_data(shot,maintracename,timebase,nverbose);
-%
-% since depends on shot number for using auto fit and thomson or thomson edge, use tracename and function here
-%
-% maintracename = 'thomson' or 'thomson_edge
+function [firthomratio] = get_fir_thom_rat_data(shot,timebase,nverbose)
 %
 % return fir_to_thomson ratio on time=timebase
 %
-% normally should use "main" thomson one built from auto fit if possible
-% take edge one if need be
-%
-% default: maintracename = "thomson"
+% Use of autofit or thomson fir ratio depends on shot number
 %
-maintracename_eff = 'thomson';
-if exist('maintracename') && ~isempty(maintracename)
-  maintracename_eff = maintracename;
-end
+% Update on May 5th 2020 (amerle):
+% Node for edge system thomson.edge:fir_thom_rat was never used so this
+% option is now removed. The edge system existed for shots 24382:49895
+% which lies entirely in the range for the autofit ratio. Note also that a
+% separate fir ratio for edge and main makes very little sense as none of
+% them cover the whole chord height.
 
 firthomratio = NaN;
-if ~exist('timebase') || isempty(timebase)
+if ~exist('timebase','var') || isempty(timebase)
   if (nverbose>=1)
     disp('need a timebase in get_fir_thom_rat_data')
   end
   return
 end
 firthomratio = NaN*ones(size(timebase));
-
-if strcmp(maintracename_eff,'thomson')
-  if shot>=23801
-    tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
-    if isempty(tracefirrat.data) || ischar(tracefirrat.data)
-      if (nverbose>=1); disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat'); end
-      tracefirrat=tdi('\results::thomson:fir_thom_rat');
-    end
-  else
+if shot>=23801
+  tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
+  if isempty(tracefirrat.data) || ischar(tracefirrat.data)
+    if (nverbose>=1); disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat'); end
     tracefirrat=tdi('\results::thomson:fir_thom_rat');
-    if isempty(tracefirrat.data) || ischar(tracefirrat.data)
-      if (nverbose>=1); disp('problem with \results::thomson:fir_thom_rat: empty'); end
-    end
-  end
-elseif strcmp(maintracename_eff,'thomson_edge')
-  if shot>=23801
-    tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
-    if isempty(tracefirrat.data) || ischar(tracefirrat.data)
-      if (nverbose>=1); disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat'); end
-      tracefirrat=tdi('\results::thomson:fir_thom_rat');
-    end
-  else
-    tracefirrat=tdi('\results::thomson_edge:fir_thom_rat');
-    if isempty(tracefirrat.data) || ischar(tracefirrat.data)
-      if (nverbose>=1); disp('problem with \results::thomson_edge:fir_thom_rat: empty'); end
-    end
   end
 else
-  if (nverbose>=1); disp('bad input in get_fir_thom_rat_data'); end
-  return
+  tracefirrat=tdi('\results::thomson:fir_thom_rat');
+  if isempty(tracefirrat.data) || ischar(tracefirrat.data)
+    if (nverbose>=1); disp('problem with \results::thomson:fir_thom_rat: empty'); end
+  end
 end
 
 if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) && any(isfinite(tracefirrat.data)) ...
@@ -3461,7 +3428,7 @@ if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) && any(isfinite(trace
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,doedge,nverbose);
+function [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,systems,nverbose)
 %
 try
   time=mdsdata('\results::thomson:times');
@@ -3481,73 +3448,162 @@ catch
   return
 end
 if isempty(time) || ischar(time)
-  thomsontimes=time;
   if (nverbose>=1) && shot<100000
     warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
     disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
   end
   return
 end
-edge_str_ = '';
-edge_str_dot = '';
-if doedge
-  edge_str_ = '_edge';
-  edge_str_dot = '.edge';
+
+systems = cellstr(systems);
+if numel(systems) == 1 && strcmpi(systems{1},'all')
+  if shot > 24381 && shot < 49896, systems = {'main','edge'};
+  else,                            systems = {'main'};
+  end
 end
+nsys = numel(systems);
+data_tmp = repmat(struct,nsys,1);
+status = false(nsys,1);
+
+for is = 1:nsys
+  switch lower(systems{is})
+    case 'main'
+      edge_str = '';
+    case 'edge'
+      edge_str = '.edge';
+    otherwise
+      warning('Unrecognised thomson scattering system name: %s',systems{is});
+      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+      return
+  end
+  
+  nodenameeff = ['\results::thomson' edge_str ':' data_request_eff(1:2)];
+  tracetdi=tdi(nodenameeff);
+  if isempty(tracetdi.data) || ischar(tracetdi.data) || isempty(tracetdi.dim)
+    data_tmp(is).data = [];
+    data_tmp(is).error_bar = [];
+    data_tmp(is).data_fullpath = '';
+    data_tmp(is).dim = {[],[]};
+    continue
+  end
+  
+  data_tmp(is).data=tracetdi.data.'; % Thomson data as (t,z)
+  tracestd=tdi(['\results::thomson'  edge_str ':' data_request_eff(1:2) ':error_bar']);
+  data_tmp(is).error_bar=tracestd.data';
+  data_tmp(is).data_fullpath=nodenameeff;
+  z=mdsdata(['\diagz::thomson_set_up' edge_str ':vertical_pos']);
+  data_tmp(is).dim={z,time};
+  data_tmp(is).units=tracetdi.units;
+  data_tmp(is).system=repmat({lower(systems{is})},numel(z),1);
+  status(is) = true;
 
-nodenameeff = ['\results::thomson' edge_str_dot ':' data_request_eff(1:2)];
-tracetdi=tdi(nodenameeff);
-if isempty(tracetdi.data) || ischar(tracetdi.data) || isempty(tracetdi.dim)
-  gdat_data.error_bar = [];
-  gdat_data.firrat = [];
-  gdat_data.data_raw = [];
-  gdat_data.error_bar_raw = [];
-  return
 end
 
-gdat_data.data=tracetdi.data'; % Thomson data as (t,z)
-tracestd=tdi(['\results::thomson'  edge_str_dot ':' data_request_eff(1:2) ':error_bar']);
-gdat_data.error_bar=tracestd.data';
-gdat_data.data_fullpath=[nodenameeff '; error_bar'];
-% add fir if ne requested
+% One system with good data
+status_any = any(status);
+isys_ref = find(status,1,'first');
+
+gdat_data.data = vertcat(data_tmp.data);
+gdat_data.error_bar = vertcat(data_tmp.error_bar);
+if status_any,gdat_data.data_fullpath = [strjoin({data_tmp(status).data_fullpath},'; '),'; error_bar'];end
+z = arrayfun(@(x) x.dim{1}, data_tmp, 'UniformOutput', false);
+gdat_data.dim = {vertcat(z{:}),time};
+gdat_data.dimunits = {'Z [m]' ; 'time [s]'};
+gdat_data.x = gdat_data.dim{1};
+gdat_data.t = time;
+gdat_data.system=vertcat(data_tmp.system);
+if status_any,gdat_data.units = data_tmp(isys_ref).units;end
+
+% add fir if ne requested 
+% NOTE: not a system dependent value so far. Only the main fir ratio was ever stored
 if strcmp(data_request_eff(1:2),'ne')
-  tracefirrat_data = get_fir_thom_rat_data(shot,['thomson' edge_str_],time,nverbose);
+  tracefirrat_data = get_fir_thom_rat_data(shot,time,nverbose);
   gdat_data.firrat=tracefirrat_data;
   gdat_data.data_raw = gdat_data.data;
-  gdat_data.data = gdat_data.data_raw * diag(tracefirrat_data);
   gdat_data.error_bar_raw = gdat_data.error_bar;
-  gdat_data.error_bar = gdat_data.error_bar_raw * diag(tracefirrat_data);
+  if status_any
+    gdat_data.data = gdat_data.data_raw * diag(tracefirrat_data);
+    gdat_data.error_bar = gdat_data.error_bar_raw * diag(tracefirrat_data);
+  end
   gdat_data.data_fullpath=[gdat_data.data_fullpath '; fir_thom_rat ; _raw without firrat'];
-  ij=find(isfinite(tracefirrat_data));
-  if isempty(ij)
+  if ~any(isfinite(tracefirrat_data))
     gdat_data.data_fullpath=[gdat_data.data_fullpath ' WARNING use ne_raw since firrat has NaNs thus ne=NaNs; fir_thom_rat ; _raw without firrat'];
     disp('***********************************************************************')
     disp('WARNING: ne from Thomson has fir_thom_rat with Nans so only ne_raw can be used');
     disp('***********************************************************************')
   end
 end
-z=mdsdata(['\diagz::thomson_set_up' edge_str_dot ':vertical_pos']);
-gdat_data.dim=[{z};{time}];
-gdat_data.dimunits=[{'Z [m]'} ; {'time [s]'}];
-gdat_data.x=z;
-gdat_data.t=time;
-% isfield does not work since tracetdi is not a 'struct' but a tdi object, thus use fieldnames
-if any(strcmp(fieldnames(tracetdi),'units'))
-  gdat_data.units=tracetdi.units;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [psiscatvol,psi_max] = get_thomson_psiscatvol(shot,systems,substr_liuqe,nverbose)
+% Subfunction to join multiple systems if needed
+%
+systems = cellstr(systems);
+if numel(systems) == 1 && strcmpi(systems{1},'all')
+  if shot > 24381 && shot < 49896, systems = {'main','edge'};
+  else,                            systems = {'main'};
+  end
+end
+nsys = numel(systems);
+[psi_max_,psiscatvol_] = deal(cell(nsys,1));
+status = false(nsys,1);
+
+validdata = @(x) ~isempty(x.data) & ~ischar(x.data) & ~isempty(x.dim);
+
+for is = 1:nsys
+  switch lower(systems{is})
+    case 'main'
+      edge_str = '';
+    case 'edge'
+      edge_str = '.edge';
+    otherwise
+      warning('Unrecognised thomson scattering system name: %s',systems{is});
+      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+      return
+  end
+  
+  psi_max_{is}   =gdat_tcv([],['\results::thomson',edge_str,':psi_max'    substr_liuqe],'nverbose',nverbose);
+  psiscatvol_{is}=gdat_tcv([],['\results::thomson',edge_str,':psiscatvol' substr_liuqe],'nverbose',nverbose);
+  
+  status(is) = validdata(psi_max_{is}) & validdata(psiscatvol_{is});
+end
+
+% Systems with good data
+nsysv = sum(status);
+
+if nsysv == 0
+  psi_max = psi_max_{1};
+  psiscatvol = psiscatvol_{1};
+else
+  % Select and concatenate
+  psi_max_v = [psi_max_{status}];
+  psiscatvol_v = [psiscatvol_{status}];
+  
+  psi_max = psi_max_v(1);
+  psi_max.data = horzcat(psi_max_v.data);
+  psi_max.data_fullpath = strjoin({psi_max_v.data_fullpath},'; ');
+  psi_max.label = strjoin({psi_max_v.data_fullpath},'; ');
+  
+  psiscatvol = psiscatvol_v(1);
+  psiscatvol.data = horzcat(psiscatvol_v.data);
+  z = arrayfun(@(x) x.dim{2}, psiscatvol_v, 'UniformOutput', false);
+  psiscatvol.dim{2} = vertcat(z{:});
+  psiscatvol.t = psiscatvol.dim{2}; %TODO: Time and space are inverted for now
+  psiscatvol.data_fullpath = strjoin({psiscatvol_v.data_fullpath},'; ');
+  psiscatvol.label = strjoin({psiscatvol_v.data_fullpath},'; ');
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, system_str_dot, psitbx_str, gdat_params )
+function [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, psitbx_str, gdat_params )
 % Input arguments:
 %   gdat_data: current structure containing TS data
 %   zshift: vertical shift of the equilibrium, can be a vector for time-dependent data
-%   system_str_dot: '' or '.edge' indicating thomson syste,
 %   psitbx_str: source for the equilibrium map, psitbx-style
 %   gdat_params: parameters for the gdat call % TODO: Why is it different from gdat_data.gdat_params?
 
 % Use psitbx
 t_th = gdat_data.t;
-th_z = mdsdata(['dim_of(\thomson' system_str_dot ':te,1)']); % TODO: Isn't this gdat_data.dim{1}?
+th_z = gdat_data.x;
 th_r = 0.9*ones(size(th_z));
 ntt = numel(t_th);
 nch = numel(th_z); % number of chords