diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index b9ae1b6f5655f1938d46419175b47b6811d08148..17fc6996771adc33d57106bfa80a8d9dc7588db4 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/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,33 +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_dot = '';
- if gdat_data.gdat_params.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, psitbx_str, gdat_params);
-
- % Add the results to the output of gdat
- gdat_data.psiscatvol = psiscatvol;
- gdat_data.psi_max = psi_max;
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;
@@ -2018,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 ...
@@ -2075,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
@@ -2089,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
@@ -2504,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
@@ -2633,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 ''')'];
@@ -3435,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');
@@ -3455,57 +3448,149 @@ 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 = '';
-if doedge
- edge_str = '.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 ':' 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 ':' 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,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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%