diff --git a/crpptbx_new/TCV/gdat_tcv.m b/crpptbx_new/TCV/gdat_tcv.m
index 01ab654d6b386a43b01029d737c9c6f5415040de..536caf2293128e58cdd69b25e55f587b1dd3fa86 100644
--- a/crpptbx_new/TCV/gdat_tcv.m
+++ b/crpptbx_new/TCV/gdat_tcv.m
@@ -299,20 +299,20 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi')
gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
gdat_data.units = aatmp.units;
gdat_data.dimunits = aatmp.dimunits;
- if mapping_for_tcv.new_timedim>0 && mapping_for_tcv.new_timedim ~= mapping_for_tcv.timedim
- % shift timedim to new_timedim data(i,j,...itime,k,...) -> data(i,inewtime,j,...,k,...)
+ if mapping_for_tcv.gdat_timedim>0 && mapping_for_tcv.gdat_timedim ~= mapping_for_tcv.timedim
+ % shift timedim to gdat_timedim data(i,j,...itime,k,...) -> data(i,inewtime,j,...,k,...)
% note that this means that gdat_data.x and gdat_data.t are same and correct,
% only .data, .dim and .dimunits need to be changed
iprev=[1:nbdims];
- ij=find(dim_nontim>mapping_for_tcv.new_timedim-1);
- inew=[1:mapping_for_tcv.new_timedim-1 mapping_for_tcv.timedim dim_nontim(ij)];
+ ij=find(dim_nontim>mapping_for_tcv.gdat_timedim-1);
+ inew=[1:mapping_for_tcv.gdat_timedim-1 mapping_for_tcv.timedim dim_nontim(ij)];
data_sizes = size(aatmp.data);
gdat_data.data = NaN*ones(data_sizes(inew));
abcol=ones(1,nbdims)*double(':'); abcomma=ones(1,nbdims)*double(',');
dimstr_prev=['(' repmat(':,',1,mapping_for_tcv.timedim-1) 'it,' ...
repmat(':,',1,nbdims-mapping_for_tcv.timedim-1) ':)'];
- dimstr_new=['(' repmat(':,',1,mapping_for_tcv.new_timedim-1) 'it,' ...
- repmat(':,',1,nbdims-mapping_for_tcv.new_timedim-1) ':)'];
+ dimstr_new=['(' repmat(':,',1,mapping_for_tcv.gdat_timedim-1) 'it,' ...
+ repmat(':,',1,nbdims-mapping_for_tcv.gdat_timedim-1) ':)'];
% eval gdat_data.data(;,:,...,it,...) = aatmp.data(:,:,:,it,...);
for it=1:size(aatmp.data,mapping_for_tcv.timedim)
shift_eval = ['gdat_data.data' dimstr_new ' = aatmp.data' dimstr_prev ';'];
@@ -320,6 +320,8 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi')
end
gdat_data.dim = aatmp.dim(inew);
gdat_data.dimunits = aatmp.dimunits(inew);
+ else
+ mapping_for_tcv.gdat_timedim = mapping_for_tcv.timedim;
end
gdat_data.data_fullpath=[mapping_for_tcv.expression substr_tdi];
@@ -421,10 +423,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
gdat_data.units = tracetdi.units;
end
gdat_data.dimunits = tracetdi.dimunits;
- gdat_data.request_description = ['vacuum magnetic field at R0=' num2str(R0EXP) 'm'];
+ gdat_data.request_description = ['vacuum magnetic field at R0=' num2str(R0EXP) 'm; COCOS=17'];
case {'betan'}
- % 100*beta / Ip[MA] * B0[T] * a[m]
+ % 100*beta / |Ip[MA] * B0[T]| * a[m]
% get B0 from gdat_tcv, without re-opening the shot and using the same parameters except data_request
% easily done thanks to structure call for options
params_eff = gdat_data.gdat_params;
@@ -450,14 +452,90 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
b0_t = interp1(b0.dim{1}(ij),b0.data(ij),gdat_data.t);
ij=find(~isnan(a_minor.data));
a_minor_t = interp1(a_minor.dim{1}(ij),a_minor.data(ij),gdat_data.t);
- gdat_data.data = 100.*beta.data ./ ip_t.*1.e6 .* b0_t .* a_minor_t;
+ gdat_data.data = 100.*beta.data ./ abs(ip_t).*1.e6 .* abs(b0_t) .* a_minor_t;
+ gdat_data.data_fullpath='100*beta/ip*1e6*b0*a_minor, each from gdat_tcv';
+ gdat_data.units = '';
+ gdat_data.dimunits = beta.dimunits;
+
+ case {'cxrs'}
+ %not yet finished, just started
+ return
+ % load typical data from cxrs, Ti, ni, vtori and vpoli (if available), as well as zeff from cxrs
+ % if 'fit' option is added: 'fit',1, then the fitted profiles are returned
+ %
+ sub_nodes = {'Ti','vi_tor','vi_pol','ni','zeff'}; % first node is also copied into data, choose "default' one
+ % sub_nodes_fit = {'Tifit','vi_torfit','vi_polfit','nifit','zefffit'};
+ params_eff = gdat_data.gdat_params;
+ % use A. Karpushov routine to get profiles and then copy the data or the fitted profiles
+ param_cxrs.k_plot=0; param_cxrs.k_debug=0;
+ cxrs_profiles = CXRS_get_profiles(48836,[],[],param_cxrs);
+ if isfield(params_eff,'fit') && params_eff.fit>0
+ sub_nodes_eff = sub_nodes_fit;
+ else
+ sub_nodes_eff = sub_nodes;
+ end
+
+
+ gdat_data.dim = beta.dim;
+ gdat_data.t = beta.dim{1};
+ gdat_data.data = beta.data;
+ ij=find(~isnan(ip.data));
+ ip_t = interp1(ip.dim{1}(ij),ip.data(ij),gdat_data.t);
+ ij=find(~isnan(b0.data));
+ b0_t = interp1(b0.dim{1}(ij),b0.data(ij),gdat_data.t);
+ ij=find(~isnan(a_minor.data));
+ a_minor_t = interp1(a_minor.dim{1}(ij),a_minor.data(ij),gdat_data.t);
+ gdat_data.data = 100.*beta.data ./ abs(ip_t).*1.e6 .* abs(b0_t) .* a_minor_t;
gdat_data.data_fullpath='100*beta/ip*1e6*b0*a_minor, each from gdat_tcv';
gdat_data.units = '';
gdat_data.dimunits = beta.dimunits;
+ case {'eqdsk'}
+ %
+ time=1.; % default time
+ if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time)
+ time = gdat_data.gdat_params.time;
+ end
+ gdat_data.gdat_params.time = time;
+ gdat_data.t = time;
+ zshift = 0.;
+ if isfield(gdat_data.gdat_params,'zshift') && ~isempty(gdat_data.gdat_params.zshift)
+ zshift = gdat_data.gdat_params.zshift;
+ end
+ gdat_data.gdat_params.zshift = zshift;
+ % use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2
+ [fnames_readresults]=read_results_for_chease(shot,time,liuqe_version,3,[],[],[],zshift,0);
+ eqdskval=read_eqdsk(fnames_readresults{4},7); % LIUQE is 17 but read_results divided psi by 2pi thus 7
+ for i=1:length(fnames_readresults)
+ unix(['rm ' fnames_readresults{i}]);
+ end
+ % transform to cocos=2 since read_results originally assumed it was cocos=2
+ cocos_in = 2;
+ [eqdsk_cocos_in, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdskval,[7 cocos_in]);
+ fnamefull = fullfile(['/tmp/' getenv('USER')],['EQDSK_' num2str(shot) 't' num2str(time,'%.4f')]);
+ % We still write COCOS=2 case, since closer to standard (in /tmp)
+ write_eqdsk(fnamefull,eqdsk_cocos_in,cocos_in);
+ % Now gdat_tcv should return the convention from LIUQE which is COCOS=17, except if specified in option
+ % create standard filename name from shot, time (cocos will be added by write_eqdsk)
+ cocos_out = 17;
+ if isfield(gdat_data.gdat_params,'cocos') && ~isempty(gdat_data.gdat_params.cocos)
+ cocos_out = gdat_data.gdat_params.cocos;
+ end
+ [eqdsk_cocosout, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdsk_cocos_in,[cocos_in cocos_out]);
+ gdat_data.eqdsk = write_eqdsk(fnamefull,eqdsk_cocosout,cocos_out);
+ gdat_data.data = gdat_data.eqdsk.psi;
+ gdat_data.dim{1} = gdat_data.eqdsk.rmesh;
+ gdat_data.dim{2} = gdat_data.eqdsk.zmesh;
+ gdat_data.dim{3} = gdat_data.t;
+ gdat_data.x = gdat_data.dim(1:2);
+ gdat_data.data_fullpath=['psi(R,Z) and eqdsk from read_eqdsk from LIUQE' num2str(liuqe_version) ';zshift=' num2str(zshift)];
+ gdat_data.units = 'T m^2';
+ gdat_data.dimunits = {'m','m','s'};
+ gdat_data.request_description = ['data=psi, x=(R,Z), eqdsk contains eqdsk structure with which ' ...
+ 'plot_eqdsk, write_eqdsk, read_eqdsk can be used'];
+
case {'ne','te'}
% ne or Te from Thomson data on raw z mesh vs (z,t)
- mdsopen(shot);
nodenameeff=['\results::thomson:' data_request_eff];
tracetdi=tdi(nodenameeff);
tracestd=tdi(['\results::thomson:' data_request_eff ':error_bar']);
@@ -480,7 +558,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
return
end
if strcmp(data_request_eff(1:2),'ne')
- tracefirrat_data = get_fir_thom_rat_data('thomson',time);
+ tracefirrat_data = get_fir_thom_rat_data(shot,'thomson',time);
gdat_data.data_abs = gdat_data.data * diag(tracefirrat_data);
gdat_data.error_bar_abs = gdat_data.error_bar * diag(tracefirrat_data);
gdat_data.firrat=tracefirrat_data;
@@ -513,13 +591,13 @@ end
if ishot==shot; mdsclose; end
-gdat_data.mapping_for_tcv = mapping_for_tcv;
+gdat_data.mapping_for.tcv = mapping_for_tcv;
error_status=0;
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function [firthomratio] = get_fir_thom_rat_data(maintracename,timebase);
+function [firthomratio] = get_fir_thom_rat_data(shot,maintracename,timebase);
%
% since depends on shot number for using auto fit and thomson or thomson edge, use tracename and function here
%
diff --git a/crpptbx_new/TCV/tcv_requests_mapping.m b/crpptbx_new/TCV/tcv_requests_mapping.m
index 435b02ae4a27479304ed5a83161d0fcd791ddb31..465d3ab15cc42182542a74834eca57a81cc0ea92 100644
--- a/crpptbx_new/TCV/tcv_requests_mapping.m
+++ b/crpptbx_new/TCV/tcv_requests_mapping.m
@@ -5,10 +5,12 @@ mapping = struct(...
'label', '', ...
'method', '', ...
'expression','', ...
- 'timedim', -1, ... % dim which is the time, to copy in .t, the other dims are in .x (-1 means last dimension)
- 'new_timedim',0, ... % if need to reshape data and dim orders to have timedim as new_timedim (shifting time to new_timedim)
+ 'timedim', -1, ... % dim which is the time is the database, to copy in .t, the other dims are in .x (-1 means last dimension)
+ 'gdat_timedim',[], ... % if need to reshape data and dim orders to have timedim as gdat_timedim (shifting time to gdat_timedim)
'min', -inf, ...
'max', inf);
+% Note that gdat_timedim is set to timedim at end of this function if empty
+% gdat_timedim should have effective index of the time dimension in gdat
if ~exist('data_request') || isempty(data_request)
return
@@ -27,140 +29,178 @@ mapping.label = data_request;
% label is used for plotting
switch lower(data_request)
case 'a_minor'
+ mapping.timedim = 1;
mapping.label = 'a(LCFS)';
mapping.method = 'switchcase';
mapping.expression = '';
case 'b0'
+ mapping.timedim = 1;
mapping.label = 'B_0';
mapping.method = 'switchcase';
mapping.expression = '';
case 'beta'
+ mapping.timedim = 1;
mapping.label = '\beta';
mapping.method = 'tdiliuqe';
mapping.expression = '\results::beta_tor';
case 'betan'
+ mapping.timedim = 1;
mapping.label = '\beta_N';
mapping.method = 'switchcase';
mapping.expression = '';
case 'betap'
+ mapping.timedim = 1;
mapping.label = '\beta_p';
mapping.method = 'tdiliuqe';
mapping.expression = '\results::beta_pol';
case 'cxrs'
+ mapping.timedim = 2;
mapping.label = 'cxrs';
mapping.method = 'switchcase';
mapping.expression = '';
case 'delta'
+ mapping.timedim = 1;
mapping.method = 'tdiliuqe';
mapping.expression = '\results::delta_edge';
% mapping.method = 'function';
% mapping.expression = ['tdi(''\results::q_psi'');']; % for tests
case 'delta_top'
+ mapping.timedim = 1;
mapping.label = 'delta\_top';
mapping.method = 'tdiliuqe';
mapping.expression = '\results::delta_ed_top';
case 'delta_bottom'
+ mapping.timedim = 1;
mapping.label = 'delta\_bottom';
mapping.method = 'tdiliuqe';
mapping.expression = '\results::delta_ed_bot';
case 'ece'
+ mapping.timedim = 2;
mapping.method = 'switchcase';
mapping.expression = '';
case 'eqdsk'
+ mapping.timedim = 0;
mapping.method = 'switchcase'; % could use function make_eqdsk directly?
mapping.expression = '';
case 'halpha'
+ mapping.timedim = 1;
mapping.label = 'Halpha';
mapping.method = 'tdi';
mapping.expression = '\base::pd:pd_011';
case 'ioh'
+ mapping.timedim = 1;
mapping.label = 'I ohmic transformer';
mapping.method = 'tdi';
mapping.expression = [{'\magnetics::ipol[*,$1]'} {'OH_001'}];
case 'ip'
+ mapping.timedim = 1;
mapping.label = 'Plasma current';
mapping.method = 'tdi';
mapping.expression = '\magnetics::iplasma:trapeze';
case 'kappa'
+ mapping.timedim = 1;
mapping.method = 'tdiliuqe';
mapping.expression = '\results::kappa_edge';
case 'mhd'
+ mapping.timedim = 1;
mapping.label = 'n=1,2, etc';
mapping.method = 'switchcase';
mapping.expression = '';
case 'ne'
+ mapping.timedim = 2;
mapping.method = 'switchcase';
case 'neint'
+ mapping.timedim = 1;
mapping.label = 'line integrated el. density';
mapping.method = 'tdi';
mapping.expression = '\results::fir:lin_int_dens';
case 'nel'
+ mapping.timedim = 1;
mapping.label = 'line-averaged el. density';
mapping.method = 'tdi';
mapping.expression = '\results::fir:n_average';
case 'nerho'
+ mapping.timedim = 2;
mapping.label = 'ne';
mapping.method = 'switchcase';
case 'neterho'
+ mapping.timedim = 2;
mapping.label = 'ne and Te';
mapping.method = 'switchcase';
case 'ni'
+ mapping.timedim = 2;
mapping.method = 'switchcase'; % especially since might have option fit, etc
case 'powers'
mapping.label = 'various powers';
mapping.method = 'switchcase';
case 'q0'
+ mapping.timedim = 1;
mapping.method = 'tdiliuqe';
mapping.expression = '\results::q_zero';
case 'q95'
+ mapping.timedim = 1;
mapping.method = 'tdiliuqe';
mapping.expression = '\results::q_95';
case 'qedge'
+ mapping.timedim = 1;
mapping.method = 'tdiliuqe';
mapping.expression = '\results::q_edge';
case 'qrho'
+ mapping.timedim = 2;
mapping.label = 'q';
mapping.method = 'switchcase';
case 'rgeom'
+ mapping.timedim = 1;
mapping.label = 'Rgeom';
mapping.method = 'switchcase';
case 'rhovol'
+ mapping.timedim = 2;
mapping.label = 'rhovol\_norm';
mapping.method = 'switchcase'; % from conf if exist otherwise computes it
case 'rmag'
+ mapping.timedim = 1;
mapping.label = 'R\_magaxis';
mapping.method = 'tdiliuqe';
mapping.expression = '\results::r_axis';
case 'sxr'
+ mapping.timedim = 2;
mapping.method = 'switchcase';
case 'te'
+ mapping.timedim = 2;
mapping.label = 'Te';
mapping.method = 'switchcase';
case 'terho'
+ mapping.timedim = 2;
mapping.label = 'Te';
mapping.method = 'switchcase';
case 'ti'
+ mapping.timedim = 2;
mapping.label = 'Ti';
mapping.method = 'switchcase';
case 'transp'
mapping.label = 'transp output';
mapping.method = 'switchcase';
case 'vloop'
+ mapping.timedim = 1;
mapping.label = '';
mapping.method = 'tdi';
- mapping.expression = '';
+ mapping.expression = [{'\magnetics::vloop[*,$1]'} {'001'}];
case 'vol'
+ mapping.timedim = 2;
mapping.label = 'Volume';
mapping.method = 'switchcase';
% mapping.expression = '\results::psitbx:vol'; (if exists for liuqe2 and 3 as well)
case 'zeff'
+ mapping.timedim = 1;
mapping.label = 'zeff from Ip-Ibs';
mapping.method = 'tdi';
mapping.expression = '\results::ibs:z_eff';
case 'zgeom'
+ mapping.timedim = 1;
mapping.label = 'Zgeom';
mapping.method = 'switchcase';
case 'zmag'
+ mapping.timedim = 1;
mapping.label = 'Zmagaxis';
mapping.method = 'tdiliuqe';
mapping.expression = '\results::z_axis';
@@ -175,4 +215,6 @@ switch lower(data_request)
end
-
+if isempty(mapping.gdat_timedim)
+ mapping.gdat_timedim = mapping.timedim;
+end
diff --git a/crpptbx_new/gdat_plot.m b/crpptbx_new/gdat_plot.m
index 02c616e88bd508ecd80afec3aea57fb760650944..04a3c9b4c4a07e91d50ea8400cd1ccfc36019fb7 100644
--- a/crpptbx_new/gdat_plot.m
+++ b/crpptbx_new/gdat_plot.m
@@ -25,10 +25,17 @@ if prod(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty
fighandle = figure(abs(gdat_data.gdat_params.doplot));
hold all
end
- plot(gdat_data.t,gdat_data.data);
- title([gdat_data.gdat_params.machine ' #' num2str(gdat_data.shot)]);
- ylabel(gdat_data.label);
- zoom on;
+ if any(find(size(gdat_data.data)==length(gdat_data.t)))
+ plot(gdat_data.t,gdat_data.data);
+ title([gdat_data.gdat_params.machine ' #' num2str(gdat_data.shot)]);
+ if isfield(gdat_data,'mapping_for')
+ xlabel(['time [' gdat_data.dimunits{gdat_data.mapping_for.(gdat_data.gdat_params.machine).gdat_timedim} ']']);
+ else
+ xlabel(['time']);
+ end
+ ylabel([gdat_data.label '[' gdat_data.units ']']);
+ zoom on;
+ end
else
disp('cannot plot gdat_data, has empty data or t field')
end