From f2a0b985d29dadbcee812b5e814bccdeb36a30bc Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Sun, 27 Sep 2015 16:20:15 +0000
Subject: [PATCH] added equil, eqdsk, ne, te, ne_rho, te_rho and direct
aug_mapping cases in new gdat for AUG, improved geteqdskAUG as well
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@5040 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/AUG/aug_requests_mapping.m | 247 ++++++++----
crpptbx/AUG/gdat_aug.m | 613 ++++++++++++++++++++++++-----
crpptbx/AUG/geteqdskAUG.m | 54 ++-
crpptbx/AUG/loadAUGdata.m | 2 +-
crpptbx/AUG/rdaAUG_eff.m | 3 +-
crpptbx/TCV/gdat_tcv.m | 20 +-
crpptbx/TCV/tcv_requests_mapping.m | 1 +
crpptbx/gdat_plot.m | 23 +-
8 files changed, 753 insertions(+), 210 deletions(-)
diff --git a/crpptbx/AUG/aug_requests_mapping.m b/crpptbx/AUG/aug_requests_mapping.m
index 15b98dd6..17dc0dff 100644
--- a/crpptbx/AUG/aug_requests_mapping.m
+++ b/crpptbx/AUG/aug_requests_mapping.m
@@ -1,14 +1,19 @@
function mapping = aug_requests_mapping(data_request)
+%
+% Information pre-defined for gdat_aug, you can add cases here to match official cases in gdat_aug, allowing backward compatibility
+%
% Defaults
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
@@ -18,78 +23,140 @@ end
mapping.label = data_request;
% for AUG, following choices are set so far:
-% method = 'tdi' and then expression is the string within tdi (usual case when there is a direct link to an existing signal)
-% with tdi, if expression cell array, call tdi(cell{1},cell{2},...)
-% method = 'function', then expression is the funtion to call which should return the correct structure
+% method = 'signal' then expression contains the shotfile, diagnostic and if needed the experiment
+% expression is a cell array
+% method = 'expression', then expression is the expression to return gdat_tmp...
% method = 'switchcase', then there will be a specific case within gdat_aug (usual case when not directly a signal)
%
% label is used for plotting
+if iscell(data_request) % || (~ischar(data_request) && length(data_request)>1)
+ mapping.label = data_request;
+ mapping.method = 'signal'; % assume a full tracename is given, so just try with tdi (could check there are some ":", etc...)
+ mapping.expression = data_request;
+ mapping.gdat_timedim = mapping.timedim;
+ return
+end
+
switch lower(data_request)
+ case 'a_minor'
+ mapping.timedim = 1;
+ mapping.label = 'a\_minor';
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''r_inboard'';' ...
+ 'gdat_tmp=gdat_aug(shot,params_eff);gdat_tmp.r_inboard=gdat_tmp.data;' ...
+ 'params_eff.data_request=''r_outboard'';' ...
+ 'gdat_tmp2=gdat_aug(shot,params_eff);gdat_tmp.r_outboard=gdat_tmp2.data;' ...
+ 'gdat_tmp.data = 0.5.*(gdat_tmp2.data-gdat_tmp.data);gdat_tmp.label=''' mapping.label ''';' ...
+ 'gdat_tmp.gdat_request=''' data_request ''';'];
case 'b0'
+ mapping.timedim = 1;
mapping.label = 'B_0';
- mapping.method = 'switchcase';
- mapping.expression = '';
+ mapping.method = 'signal';
+ mapping.expression = [{'FPC'},{'BTF'}];
+ case 'beta'
+ mapping.timedim = 1;
+ mapping.label = '\beta';
+ mapping.method = 'signal';
+ mapping.expression = [{'TOT'},{'beta'}];
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''betan'';' ...
+ 'gdat_tmp=gdat_aug(shot,params_eff);' ...
+ 'params_eff.data_request=''ip'';gdat_tmp2=gdat_aug(shot,params_eff);' ...
+ 'params_eff.data_request=''b0'';gdat_tmp3=gdat_aug(shot,params_eff);' ...
+ 'params_eff.data_request=''a_minor'';gdat_tmp4=gdat_aug(shot,params_eff);' ...
+ 'tmp_data_ip=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ...
+ 'tmp_data_b0=interp1(gdat_tmp3.t,gdat_tmp3.data,gdat_tmp.t,[],NaN);' ...
+ 'tmp_data_a=interp1(gdat_tmp4.t,gdat_tmp4.data,gdat_tmp.t,[],NaN);' ...
+ 'gdat_tmp.data = 0.01.*abs(gdat_tmp.data.*tmp_data_ip./1e6./tmp_data_a./tmp_data_b0);' ...
+ 'gdat_tmp.label=''' mapping.label ''';gdat_tmp.gdat_request=''' data_request ''';'];
case 'betan'
+ mapping.timedim = 1;
mapping.label = '\beta_N';
- mapping.method = 'switchcase';
- mapping.expression = '';
+ mapping.method = 'signal';
+ mapping.expression = [{'TOT'},{'beta_N'}];
case 'betap'
+ mapping.timedim = 1;
mapping.label = '\beta_p';
- mapping.method = 'tdi';
- mapping.expression = '\results::beta_pol';
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'betpol'}];
case 'cxrs'
+ mapping.timedim = 2;
mapping.label = 'cxrs';
mapping.method = 'switchcase';
mapping.expression = '';
case 'delta'
- mapping.method = 'tdi';
- mapping.expression = '\results::delta_edge';
- mapping.method = 'function';
- mapping.expression = ['tdi(''\results::q_psi'');'];
+ mapping.timedim = 1;
+ mapping.label = 'delta';
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''delta_bottom''; ' ...
+ 'gdat_tmp=gdat_aug(shot,params_eff);params_eff.data_request=''delta_top'';' ...
+ 'gdat_tmp2=gdat_aug(shot,params_eff);gdat_tmp.data = 0.5.*(gdat_tmp.data+gdat_tmp2.data);'];
case 'delta_top'
mapping.label = 'delta\_top';
- mapping.method = 'tdi';
- mapping.expression = '\results::delta_ed_top';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'delRoben'}];
case 'delta_bottom'
mapping.label = 'delta\_bottom';
- mapping.method = 'tdi';
- mapping.expression = '\results::delta_ed_bot';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'delRuntn'}];
case 'ece'
+ mapping.timedim = 2;
mapping.method = 'switchcase';
mapping.expression = '';
case 'eqdsk'
+ mapping.timedim = 2;
+ mapping.method = 'switchcase'; % could use function make_eqdsk directly?
+ mapping.expression = '';
+ case 'equil'
+ mapping.timedim = 3;
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';
+ mapping.method = 'signal';
+ mapping.expression = [{'POT'},{'ELMa-Han'}];
case 'ioh'
+ mapping.timedim = 1;
mapping.label = 'I ohmic transformer';
- mapping.method = 'tdi';
- mapping.expression = [{'\magnetics::ipol[*,$1]'} {'OH_001'}];
+ mapping.method = 'signal';
+ mapping.expression = [{'MBI'},{'IOH'}];
case 'ip'
+ mapping.timedim = 1;
mapping.label = 'Plasma current';
- mapping.method = 'tdi';
- mapping.expression = '\magnetics::iplasma:trapeze';
+ mapping.method = 'signal';
+ mapping.expression = [{'MAG'},{'Ipa'}];
case 'kappa'
- mapping.method = 'tdi';
- mapping.expression = '\results::kappa_edge';
+ mapping.timedim = 1;
+ mapping.label = '\kappa';
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'k'}];
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';
+ mapping.method = 'signal';
+ mapping.expression = [{'DCN'},{'H-1'},{'AUGD'}];
case 'nel'
+ mapping.timedim = 1;
mapping.label = 'line-averaged el. density';
- mapping.method = 'tdi';
- mapping.expression = '\results::fir:n_average';
- case 'nerho'
+ mapping.expression = [{'FPG'},{'lenH-1'},{'AUGD'}];
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''neint'';' ...
+ 'gdat_tmp=gdat_aug(shot,params_eff);params_eff.data_request=[{''FPG''},{''lenH-1''},{''AUGD''}];' ...
+ 'gdat_tmp2=gdat_aug(shot,params_eff);ij=find(gdat_tmp2.data==0);gdat_tmp2.data(ij)=NaN;' ...
+ 'tmp_data=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ...
+ 'gdat_tmp.data = gdat_tmp.data./(tmp_data+1e-5);'];
+ case 'ne_rho'
mapping.label = 'ne';
mapping.method = 'switchcase';
case 'neterho'
@@ -101,69 +168,117 @@ switch lower(data_request)
mapping.label = 'various powers';
mapping.method = 'switchcase';
case 'q0'
- mapping.method = 'tdi';
- mapping.expression = '\results::q_zero';
+ mapping.timedim = 1;
+ mapping.label = 'q_0';
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'q0'},{'AUGD'}];
case 'q95'
- mapping.method = 'tdi';
- mapping.expression = '\results::q_95';
- case 'qedge'
- mapping.method = 'tdi';
- mapping.expression = '\results::q_edge';
- case 'qrho'
+ mapping.timedim = 1;
+ mapping.label = 'q_{95}';
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'q95'},{'AUGD'}];
+ case 'q_edge'
+ mapping.timedim = 1;
+ mapping.label = 'q_{edge}}';
+ mapping.method = 'expression';
+ mapping.method = 'switchcase';
+ mapping.expression = [{'FPG'},{'q95'},{'AUGD'}];
+ mapping.expression = [];
+ case 'q_rho'
mapping.label = 'q';
mapping.method = 'switchcase';
case 'rgeom'
mapping.label = 'Rgeom';
- mapping.method = 'switchcase';
+ mapping.timedim = 1;
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''r_inboard'';' ...
+ 'gdat_tmp=gdat_aug(shot,params_eff);gdat_tmp.r_inboard=gdat_tmp.data;' ...
+ 'params_eff.data_request=''r_outboard'';' ...
+ 'gdat_tmp2=gdat_aug(shot,params_eff);gdat_tmp.r_outboard=gdat_tmp2.data;' ...
+ 'gdat_tmp.data = 0.5.*(gdat_tmp2.data+gdat_tmp.data);gdat_tmp.label=''' mapping.label ''';' ...
+ 'gdat_tmp.gdat_request=''' data_request ''';'];
+ case 'r_inboard'
+ mapping.label = 'R\_inboard';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'Rin'},{'AUGD'}];
+ case 'r_outboard'
+ mapping.label = 'R\_outboard';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'Raus'},{'AUGD'}];
case 'rhovol'
mapping.label = 'rhovol\_norm';
mapping.method = 'switchcase'; % from conf if exist otherwise computes it
case 'rmag'
mapping.label = 'R\_magaxis';
- mapping.method = 'tdi';
- mapping.expression = '\results::r_axis';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'Rmag'},{'AUGD'}];
case 'sxr'
mapping.method = 'switchcase';
case 'te'
mapping.label = 'Te';
mapping.method = 'switchcase';
- case 'terho'
+ case 'te_rho'
mapping.label = 'Te';
mapping.method = 'switchcase';
case 'ti'
mapping.label = 'Ti';
mapping.method = 'switchcase';
- case 'transp'
- mapping.label = 'transp output';
- mapping.method = 'switchcase';
case 'vloop'
- mapping.label = '';
- mapping.method = 'tdi';
- mapping.expression = '';
- case 'vol'
+ mapping.label = 'Vloop';
+ mapping.timedim = 1;
+ % mapping.method = 'signal';
+ % mapping.expression = [{'MAG'},{'ULid12'},{'AUGD'}];
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=[{''MAG''},{''ULid12''},{''AUGD''}];' ...
+ 'gdat_tmp=gdat_aug(shot,params_eff);ij=find(~isnan(gdat_tmp.data));' ...
+ 'tmp_data=interpos(gdat_tmp.t,gdat_tmp.data,-3e4);' ...
+ 'gdat_tmp.data_smooth = tmp_data;gdat_tmp.gdat_request=''vloop'';gdat_tmp.gdat_params.data_request=gdat_tmp.gdat_request;'];
+ case 'volume'
mapping.label = 'Volume';
- mapping.method = 'switchcase';
- % mapping.expression = '\results::psitbx:vol'; (if exists for liuqe2 and 3 as well)
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'Vol'},{'AUGD'}];
case 'zeff'
- mapping.label = 'zeff from Ip-Ibs';
- mapping.method = 'tdi';
- mapping.expression = '\results::ibs:z_eff';
+ mapping.label = 'zeff from cxrs';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'ZES'},{'Zeff'},{'AUGD'}];
case 'zgeom'
mapping.label = 'Zgeom';
- mapping.method = 'switchcase';
+ mapping.timedim = 1;
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=[{''FPG''},{''Zoben''},{''AUGD''}];' ...
+ 'gdat_tmp=gdat_aug(shot,params_eff);gdat_tmp.z_top=gdat_tmp.data;' ...
+ 'params_eff.data_request=[{''FPG''},{''Zunt''},{''AUGD''}];' ...
+ 'gdat_tmp2=gdat_aug(shot,params_eff);gdat_tmp.z_bottom=gdat_tmp2.data;' ...
+ 'gdat_tmp.data = 0.5.*(gdat_tmp2.data+gdat_tmp.data);gdat_tmp.label=''' mapping.label ''';' ...
+ 'gdat_tmp.gdat_request=''' data_request ''';'];
+
case 'zmag'
- mapping.label = 'Zmagaxis';
- mapping.method = 'tdi';
- mapping.expression = '\results::z_axis';
-% $$$ case ''
-% $$$ mapping.label = '';
-% $$$ mapping.method = 'tdi';
-% $$$ mapping.expression = '';
+ mapping.label = 'Z\_magaxis';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'Zmag'},{'AUGD'}];
+ %
+ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % extra AUG cases (not necessarily in official data_request name list)
+ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ case 'transp'
+ mapping.label = 'transp output';
+ mapping.method = 'switchcase';
+
+
otherwise
mapping.label = data_request;
- mapping.method = 'tdi'; % assume a full tracename is given, so just try with tdi (could check there are some ":", etc...)
+ mapping.method = 'signal'; % assume a full tracename is given, so just try with tdi (could check there are some ":", etc...)
mapping.expression = data_request;
end
-
+if isempty(mapping.gdat_timedim)
+ mapping.gdat_timedim = mapping.timedim;
+end
diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m
index f070e344..083a28ea 100644
--- a/crpptbx/AUG/gdat_aug.m
+++ b/crpptbx/AUG/gdat_aug.m
@@ -32,7 +32,7 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_aug(shot,data_req
% gdat_data.shot: shot number
% gdat_data.machine: machine providing the data
% gdat_data.gdat_request: keyword for gdat if relevant
-% gdat_data.data_fullpath: full path to the data node if known and relevant, or expression, or relevant function called if relevant
+% gdat_data.data_fullpath: full path to the data node if known and relevant, or relevant expression called if relevant
% gdat_data.gdat_params: copy gdat_params for completeness
% gdat_data.xxx: any other relevant information
%
@@ -58,6 +58,7 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_aug(shot,data_req
varargout{1}=cell(1,1);
error_status=1;
+nverbose = 1;
% construct default parameters structure
gdat_params.data_request = '';
@@ -65,9 +66,11 @@ default_machine = 'aug';
gdat_params.machine=default_machine;
gdat_params.doplot = 0;
+gdat_params.exp_name = 'AUGD';
+gdat_params.equil = 'EQI';
% construct list of keywords from global set of keywords and specific AUG set
-% get data_request names from centralized function
+% get data_request names from centralized table
data_request_names = get_data_request_names;
% add AUG specific to all:
@@ -115,8 +118,9 @@ do_mdsopen_mdsclose = 1;
if nargin>=1
if isempty(shot)
% means mdsopen(shot) already performed
- shot = mdsipmex(2,'$SHOT');
- gdat_data.shot = shot;
+% $$$ shot = mdsipmex(2,'$SHOT');
+% $$$ gdat_data.shot = shot;
+ shot=0;
do_mdsopen_mdsclose = 0;
elseif isnumeric(shot)
gdat_data.shot = shot;
@@ -144,7 +148,7 @@ if nargin>=2 && ivarargin_first_char~=1
error_status=3;
return
end
- data_request.data_request = lower(data_request.data_request);
+ %data_request.data_request = lower(data_request.data_request);
data_request_eff = data_request.data_request;
gdat_params = data_request;
else
@@ -199,8 +203,14 @@ if (nargin>=ivarargin_first_char)
end
data_request_eff = gdat_params.data_request; % in case was defined in pairs
-% if it is a request_keyword copy it:
-ij=strmatch(data_request_eff,data_request_names_all,'exact');
+% if it is a request_keyword can obtain description:
+if length(data_request_eff)==1
+ ij=strmatch(data_request_eff,data_request_names_all,'exact');
+else
+ ij=[];
+end
+gdat_data.gdat_request = data_request_eff;
+
if ~isempty(ij);
gdat_data.gdat_request = data_request_names_all{ij};
if isfield(data_request_names.all.(data_request_names_all{ij}),'description') && ~isempty(data_request_names.all.(data_request_names_all{ij}).description)
@@ -233,103 +243,115 @@ error_status = 6; % at least reached this level
mapping_for_aug = aug_requests_mapping(data_request_eff);
gdat_data.label = mapping_for_aug.label;
-ishot=NaN;
-if do_mdsopen_mdsclose
- mdsdefaultserver aug1.epfl.ch; % should be in aug general path, but set-it in the meantime...
- ishot = mdsopen(shot); % if ishot equal to shot, then mdsclose at the end
- if ishot~=shot
- warning(['cannot open shot= ' num2str(shot)])
- return
- end
-end
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% 1st treat the simplest method: "tdi"
-if strcmp(mapping_for_aug.method,'tdi')
- % need to treat liuqe2, model, etc from options....
- if iscell(mapping_for_aug.expression)
- if length(mapping_for_aug.expression)>0
- % series of arguments for tdi given in each cell
- eval_expr = ['aatmp=tdi(''' mapping_for_aug.expression{1} ''''];
- for i=2:length(mapping_for_aug.expression)
- eval_expr = [eval_expr ',''' mapping_for_aug.expression{i} ''''];
- end
- eval_expr = [eval_expr ');'];
- eval(eval_expr);
- else
- % empty or wrong expression
- error_status=701;
- return
- end
+% 1st treat the simplest method: "signal"
+if strcmp(mapping_for_aug.method,'signal')
+ exp_location = gdat_data.gdat_params.exp_name;
+ if ~isempty(mapping_for_aug.expression) && iscell(mapping_for_aug.expression) && length(mapping_for_aug.expression)>=3
+ exp_location = mapping_for_aug.expression{3};
+ elseif ~isempty(mapping_for_aug.expression) && length(mapping_for_aug.expression)>=2
+ mapping_for_aug.expression{3} = exp_location;
else
- eval_expr = ['aatmp=tdi(''' mapping_for_aug.expression ''');'];
- eval(eval_expr);
+ error_status = 101;
+ disp(['expects at least 2 cells in expression, mapping_for_aug.expression = ' mapping_for_aug.expression]);
+ return
+ end
+ gdat_data.gdat_params.exp_name = exp_location;
+ [aatmp,error_status]=rdaAUG_eff(shot,mapping_for_aug.expression{1},mapping_for_aug.expression{2},exp_location);
+ if error_status~=0
+ if nverbose>=3; disp(['error after rdaAUG in signal with data_request_eff= ' data_request_eff]); end
+ return
end
gdat_data.data = aatmp.data;
- gdat_data.dim = aatmp.dim;
- nbdims = length(gdat_data.dim);
- if mapping_for_aug.timedim==-1;
- mapping_for_aug.timedim = nbdims;
- if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_aug.timedim = nbdims-1; end
+ gdat_data.t = aatmp.t;
+ if isempty(aatmp.data)
+ return
end
- dim_nontim = setdiff([1:nbdims],mapping_for_aug.timedim);
- if ~isempty(dim_nontim)
- % since most cases have at most 2d, copy as array if data is 2D and as cell if 3D or more
- if length(dim_nontim)==1
- gdat_data.x = gdat_data.dim{dim_nontim(1)};
+ if mapping_for_aug.timedim<=0
+ % need to guess timedim
+ if prod(size(aatmp.data)) == length(aatmp.data)
+ mapping_for_aug.timedim = 1;
+ elseif length(size(aatmp.data))==2
+ % 2 true dimensions
+ if length(aatmp.t) == size(aatmp.data,1)
+ mapping_for_aug.timedim = 1;
+ else
+ mapping_for_aug.timedim = 2;
+ end
else
- gdat_data.x = gdat_data.dim(dim_nontim);
+ % more than 2 dimensions, find the one with same length to define timedim
+ mapping_for_aug.timedim = find(size(aatmp.data)==length(aatmp.t));
+ if length(mapping_for_aug.timedim); mapping_for_aug.timedim = mapping_for_aug.timedim(1); end
end
+ mapping_for_aug.gdat_timedim = mapping_for_aug.timedim
end
- gdat_data.t = gdat_data.dim{mapping_for_aug.timedim};
+ if length(size(aatmp.data))==1 | (length(size(aatmp.data))==2 & size(aatmp.data,2)==1)
+ gdat_data.dim=[{aatmp.t}];
+ gdat_data.dimunits={'time [s]'};
+ elseif length(size(aatmp.data))==2
+ gdat_data.dim{mapping_for_aug.timedim} = gdat_data.t;
+ gdat_data.dimunits{mapping_for_aug.timedim} = 'time [s]';
+ gdat_data.dim{setdiff([1:2],mapping_for_aug.timedim)} = gdat_data.x;
+ else
+ gdat_data.dim{mapping_for_aug.timedim} = gdat_data.t;
+ gdat_data.dimunits{mapping_for_aug.timedim} = 'time [s]';
+ nbdims = length(size(gdat_data.data));
+ for i=1:nbdims
+ if i~=mapping_for_aug.timedim
+ ieff=i;
+ if i>mapping_for_aug.timedim; ieff=i-1; end
+ gdat_data.dim{i} = gdat_data.x{ieff};
+ end
+ end
+ end
+ nbdims = length(gdat_data.dim);
gdat_data.units = aatmp.units;
- gdat_data.dimunits = aatmp.dimunits;
- if mapping_for_aug.new_timedim>0 && mapping_for_aug.new_timedim ~= mapping_for_aug.timedim
- % shift timedim to new_timedim data(i,j,...itime,k,...) -> data(i,inewtime,j,...,k,...)
+ if mapping_for_aug.gdat_timedim>0 && mapping_for_aug.gdat_timedim ~= mapping_for_aug.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_aug.new_timedim-1);
- inew=[1:mapping_for_aug.new_timedim-1 mapping_for_aug.timedim dim_nontim(ij)];
+ ij=find(dim_nontim>mapping_for_aug.gdat_timedim-1);
+ inew=[1:mapping_for_aug.gdat_timedim-1 mapping_for_aug.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_aug.timedim-1) 'it,' ...
repmat(':,',1,nbdims-mapping_for_aug.timedim-1) ':)'];
- dimstr_new=['(' repmat(':,',1,mapping_for_aug.new_timedim-1) 'it,' ...
- repmat(':,',1,nbdims-mapping_for_aug.new_timedim-1) ':)'];
+ dimstr_new=['(' repmat(':,',1,mapping_for_aug.gdat_timedim-1) 'it,' ...
+ repmat(':,',1,nbdims-mapping_for_aug.gdat_timedim-1) ':)'];
% eval gdat_data.data(;,:,...,it,...) = aatmp.data(:,:,:,it,...);
for it=1:size(aatmp.data,mapping_for_aug.timedim)
shift_eval = ['gdat_data.data' dimstr_new ' = aatmp.data' dimstr_prev ';'];
eval(shift_eval);
end
gdat_data.dim = aatmp.dim(inew);
- gdat_data.dimunits = aatmp.dimunits(inew);
+ % gdat_data.dimunits = aatmp.dimunits(inew);
end
gdat_data.data_fullpath=mapping_for_aug.expression;
- % end of method "tdi"
+ % end of method "signal"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-elseif strcmp(mapping_for_aug.method,'function')
- % 2nd: method="function"
- % assume expression contains function to call and which returns a structure
+elseif strcmp(mapping_for_aug.method,'expression')
+ % 2nd: method="expression"
+ % assume expression contains an expression to evaluate and which cretaes a local structure into variable gdat_tmp
% we copy the structure, to make sure default nodes are defined and to avoid if return is an closed object like tdi
- eval_expr = ['aatmp=' mapping_for_aug.expression ';'];
- eval(eval_expr);
- if isempty(aatmp) || (~isstruct(aatmp) & ~isobject(aatmp))
- warning(['function expression does not return a structure: ' eval_expr])
+ % eval_expr = [mapping_for_aug.expression ';'];
+ eval([mapping_for_aug.expression ';']);
+ if isempty(gdat_tmp) || (~isstruct(gdat_tmp) & ~isobject(gdat_tmp))
+ warning(['expression does not create a gdat_tmp structure: ' mapping_for_aug.expression])
error_status=801;
return
end
- tmp_fieldnames = fieldnames(aatmp);
- if sum(strcmp(tmp_fieldnames,'data'))==0 % note: cannot do isfield since aatmp might be an object
- warning(['function does not return a child name ''data'' for ' data_request_eff])
+ tmp_fieldnames = fieldnames(gdat_tmp);
+ if sum(strcmp(tmp_fieldnames,'data'))==0 % note: cannot do isfield since gdat_tmp might be an object
+ warning(['expression does not return a child name ''data'' for ' data_request_eff])
end
for i=1:length(tmp_fieldnames)
- gdat_data.(tmp_fieldnames{i}) = aatmp.(tmp_fieldnames{i});
+ gdat_data.(tmp_fieldnames{i}) = gdat_tmp.(tmp_fieldnames{i});
end
% add .t and .x in case only dim is provided
- % do not allow shifting of timedim since should be treated in the relevant function
+ % do not allow shifting of timedim since should be treated in the relevant expression
ijdim=find(strcmp(tmp_fieldnames,'dim')==1);
if ~isempty(ijdim)
nbdims = length(gdat_data.dim);
@@ -355,55 +377,432 @@ elseif strcmp(mapping_for_aug.method,'function')
end
gdat_data.data_fullpath=mapping_for_aug.expression;
end
- % end of method "function"
+ % end of method "expression"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_aug.method,'switchcase')
switch data_request_eff % not lower(...) since data_request_eff should be lower case already at this stage
- 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']);
- trace_fir_rat=tdi('\results::thomson:fir_thom_rat');
- gdat_data.data=tracetdi.data'; % Thomson data as (t,z)
- gdat_data.error_bar=tracestd.data';
- gdat_data.data_fullpath=[nodenameeff];
- % add correct dimensions
- try
- time=mdsdata('\results::thomson:times');
- catch
- warning('Problems with \results::thomson:times')
- disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+ case {'eqdsk'}
+ %
+ time_eqdsks=1.; % default time
+ if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time)
+ time_eqdsks = gdat_data.gdat_params.time;
+ else
+ gdat_data.gdat_params.time = time_eqdsks;
+ disp(['"time" is expected as an option, choose default time = ' num2str(time_eqdsks)]);
+ end
+ gdat_data.gdat_params.time = time_eqdsks;
+ gdat_data.t = time_eqdsks;
+ zshift = 0.;
+ if isfield(gdat_data.gdat_params,'zshift') && ~isempty(gdat_data.gdat_params.zshift)
+ zshift = gdat_data.gdat_params.zshift;
+ else
+ gdat_data.gdat_params.zshift = zshift;
+ end
+ gdat_data.gdat_params.zshift = zshift;
+ params_equil = gdat_data.gdat_params;
+ params_equil.data_request = 'equil';
+ [equil,params_equil,error_status] = gdat_aug(shot,params_equil);
+ if error_status>0
+ if nverbose>=3; disp(['problems with ' params_equil.data_request]); end
return
end
- if isempty(time) || ischar(time)
- thomsontimes=time
- warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check')
- disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+ gdat_data.gdat_params = params_equil;
+ gdat_data.equil = equil;
+ for itime=1:length(time_eqdsks)
+ time_eff = time_eqdsks(itime);
+ % use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2
+ [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(equil,time_eff,gdat_data.gdat_params.zshift);
+ eqdskAUG.fnamefull = fullfile(['/tmp/' getenv('USER')],['EQDSK_' num2str(shot) 't' num2str(time_eff,'%.4f')]);
+ cocos_out = equil.cocos;
+ if isfield(gdat_data.gdat_params,'cocos') && ~isempty(gdat_data.gdat_params.cocos)
+ cocos_out = gdat_data.gdat_params.cocos;
+ else
+ gdat_data.gdat_params.cocos = cocos_out;
+ end
+ if equil.cocos ~= cocos_out
+ [eqdsk_cocosout, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdskAUG,[equil.cocos cocos_out]);
+ else
+ eqdsk_cocosout = eqdskAUG;
+ end
+ % for several times, use array of structure for eqdsks,
+ % cannot use it for psi(R,Z) in .data and .x since R, Z might be different at different times,
+ % so project psi(R,Z) on Rmesh, Zmesh of 1st time
+ if length(time_eqdsks) > 1
+ gdat_data.eqdsk{itime} = write_eqdsk(eqdskAUG.fnamefull,eqdsk_cocosout);
+ if itime==1
+ gdat_data.data(:,:,itime) = gdat_data.eqdsk{itime}.psi;
+ gdat_data.dim{1} = gdat_data.eqdsk{itime}.rmesh;
+ gdat_data.dim{2} = gdat_data.eqdsk{itime}.zmesh;
+ else
+ xx=repmat(reshape(gdat_data.dim{1},length(gdat_data.dim{1}),1),1,size(gdat_data.eqdsk{itime}.psi,2));
+ yy=repmat(reshape(gdat_data.dim{2},1,length(gdat_data.dim{2})),size(gdat_data.eqdsk{itime}.psi,1),1);
+ aa = interpos2Dcartesian(gdat_data.eqdsk{itime}.rmesh,gdat_data.eqdsk{itime}.zmesh ...
+ ,gdat_data.eqdsk{itime}.psi,xx,yy,-1,-1);
+ gdat_data.data(:,:,itime) = aa;
+ end
+ else
+ gdat_data.eqdsk = write_eqdsk(eqdskAUG.fnamefull,eqdsk_cocosout);
+ gdat_data.data = gdat_data.eqdsk.psi;
+ gdat_data.dim{1} = gdat_data.eqdsk.rmesh;
+ gdat_data.dim{2} = gdat_data.eqdsk.zmesh;
+ end
+ end
+ 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 geteqdskAUG with ' gdat_data.gdat_params.equil ' ;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 {'equil'}
+ exp_name_eff = gdat_data.gdat_params.exp_name;
+ gdat_data.gdat_params.equil = upper(gdat_data.gdat_params.equil);
+ DIAG = gdat_data.gdat_params.equil; % 'EQI' by default at this stage, should become EQH?
+ M_Rmesh_par = sf2par(DIAG,shot,'M','PARMV');
+ M_Rmesh = M_Rmesh_par.value + 1; % nb of points
+ N_Zmesh_par = sf2par(DIAG,shot,'N','PARMV');
+ N_Zmesh = N_Zmesh_par.value + 1; % nb of points
+ NTIME_par = sf2par(DIAG,shot,'NTIME','PARMV');
+ NTIME = NTIME_par.value; % nb of points
+ Lpf_par = rdaAUG_eff(shot,DIAG,'Lpf',exp_name_eff);
+ % since June, nb of time points in EQ results is not consistent with NTIME and time
+ % It seems the first NTIME points are correct, so use this explicitely
+ NTIME_Lpf = length(Lpf_par.value);
+ if (NTIME < NTIME_Lpf)
+ if nverbose>=3; disp('WARNING: nb of times points smaller then equil results, use first NTIME points'); end
+ elseif (NTIME > NTIME_Lpf)
+ if nverbose >= 1
+ disp('ERROR: nb of times points LARGER then equil results')
+ disp('this is unexpected, so stop there and ask Olivier.Sauter@epfl.ch')
+ end
return
end
- if strcmp(data_request_eff(1:2),'ne')
- tracefirrat_data = get_fir_thom_rat_data('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;
- gdat_data.data_fullpath=[gdat_data.data_fullpath ' ; _abs includes *firrat'];
+ Lpf_tot = Lpf_par.value(1:NTIME); % nb of points: 100000*nb_SOL points + nb_core
+ Lpf_SOL = fix(Lpf_tot/100000);
+ Lpf1_t = mod(Lpf_tot,100000)+1; % nb of Lpf points
+ % since Lpf depends on time, need to load all first and then loop over time for easier mapping
+ [qpsi,e]=rdaAUG_eff(shot,DIAG,'Qpsi',exp_name_eff);
+ ndimrho = size(qpsi.data,2);
+ if ndimrho==NTIME_Lpf
+ % data seems to be transposed
+ ndimrho = size(qpsi.data,1);
+ itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
+ else
+ itotransposeback = 0;
end
- z=mdsdata('\diagz::thomson_set_up: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;
+ qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback);
+ ijnan=find(isnan(qpsi.value));
+ qpsi.value(ijnan)=0;
+ [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff);
+ psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
+ [phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',exp_name_eff);
+ phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
+ [Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',exp_name_eff);
+ Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
+ [Area,e]=rdaAUG_eff(shot,DIAG,'Area',exp_name_eff);
+ Area=adapt_rda(Area,NTIME,2*ndimrho,itotransposeback);
+ [Ri,e]=rdaAUG_eff(shot,DIAG,'Ri',exp_name_eff);
+ Ri=adapt_rda(Ri,NTIME,M_Rmesh,itotransposeback);
+ [Zj,e]=rdaAUG_eff(shot,DIAG,'Zj',exp_name_eff);
+ Zj=adapt_rda(Zj,NTIME,N_Zmesh,itotransposeback);
+ [PFM_tree,e]=rdaAUG_eff(shot,DIAG,'PFM',exp_name_eff);
+ PFM_tree=adaptPFM_rda(PFM_tree,M_Rmesh,N_Zmesh,NTIME);
+ [Pres,e]=rdaAUG_eff(shot,DIAG,'Pres',exp_name_eff);
+ Pres=adapt_rda(Pres,NTIME,2*ndimrho,itotransposeback);
+ [Jpol,e]=rdaAUG_eff(shot,DIAG,'Jpol',exp_name_eff);
+ Jpol=adapt_rda(Jpol,NTIME,2*ndimrho,itotransposeback);
+ [FFP,e]=rdaAUG_eff(shot,DIAG,'FFP',exp_name_eff);
+ if ~isempty(FFP.value)
+ FFP=adapt_rda(FFP,NTIME,ndimrho,itotransposeback);
+ else
+ FFP.value=NaN*ones(NTIME,max(Lpf1_t));
end
-
+ if strcmp(DIAG,'EQI') || strcmp(DIAG,'EQH')
+ [Rinv,e]=rdaAUG_eff(shot,DIAG,'Rinv',exp_name_eff);
+ Rinv=adapt_rda(Rinv,NTIME,ndimrho,itotransposeback);
+ [R2inv,e]=rdaAUG_eff(shot,DIAG,'R2inv',exp_name_eff);
+ R2inv=adapt_rda(R2inv,NTIME,ndimrho,itotransposeback);
+ [Bave,e]=rdaAUG_eff(shot,DIAG,'Bave',exp_name_eff);
+ Bave=adapt_rda(Bave,NTIME,ndimrho,itotransposeback);
+ [B2ave,e]=rdaAUG_eff(shot,DIAG,'B2ave',exp_name_eff);
+ B2ave=adapt_rda(B2ave,NTIME,ndimrho,itotransposeback);
+ [FTRA,e]=rdaAUG_eff(shot,DIAG,'FTRA',exp_name_eff);
+ FTRA=adapt_rda(FTRA,NTIME,ndimrho,itotransposeback);
+ else
+ Rinv.value=[]; R2inv.value=[]; Bave.value=[]; B2ave.value=[]; FTRA.value=[];
+ end
+ [LPFx,e]=rdaAUG_eff(shot,DIAG,'LPFx',exp_name_eff);
+ LPFx.value=LPFx.value(1:NTIME); LPFx.data=LPFx.value; LPFx.t=LPFx.t(1:NTIME);
+ [PFxx,e]=rdaAUG_eff(shot,DIAG,'PFxx',exp_name_eff);
+ PFxx=adapt_rda(PFxx,NTIME,max(LPFx.value)+1,itotransposeback);
+ [RPFx,e]=rdaAUG_eff(shot,DIAG,'RPFx',exp_name_eff);
+ RPFx=adapt_rda(RPFx,NTIME,max(LPFx.value)+1,itotransposeback);
+ [zPFx,e]=rdaAUG_eff(shot,DIAG,'zPFx',exp_name_eff);
+ zPFx=adapt_rda(zPFx,NTIME,max(LPFx.value)+1,itotransposeback);
+ % seems "LCFS" q-value is far too large, limit to some max (when diverted)
+ max_qValue = 30.; % Note could just put a NaN on LCFS value since ill-defined when diverted
+ for it=1:NTIME
+ Lpf1 = Lpf1_t(it);
+ % Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
+ % change it to (radial,time) and use only Lpf+1 points up to LCFS
+ ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially
+ % set NaNs to zeroes
+ if qpsi.value(ijok(1),1)<0
+ gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
+ else
+ gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
+ end
+ % get x values
+ gdat_data.psi(:,it)=psi_tree.value(it,Lpf1:-1:1)';
+ gdat_data.psi_axis(it)= gdat_data.psi(1,it);
+ gdat_data.psi_lcfs(it)= gdat_data.psi(end,it);
+ gdat_data.rhopolnorm(:,it) = sqrt(abs((gdat_data.psi(:,it)-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it))));
+ if strcmp(DIAG,'EQR');
+ % q value has only a few values and from center to edge, assume they are from central rhopol values on
+ % But they are every other point starting from 3rd
+ ijk=find(gdat_data.qvalue(:,it)~=0);
+ if length(ijk)>2
+ % now shots have non-zero axis values in eqr
+ rhoeff=gdat_data.rhopolnorm(ijk,it);
+ qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above
+ if ijk(1)>1
+ rhoeff = [0.; rhoeff];
+ qeff = [qeff(1) ;qeff];
+ end
+ ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it)));
+ qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
+ qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]);
+ else
+ qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
+ end
+ gdat_data.qvalue(:,it) = qfit;
+ end
+ % get rhotor values
+ gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)';
+ gdat_data.rhotornorm(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.phi(end,it)));
+ % get rhovol values
+ gdat_data.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)';
+ gdat_data.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
+ gdat_data.rhovolnorm(:,it) = sqrt(abs(gdat_data.vol(:,it) ./ gdat_data.vol(end,it)));
+ gdat_data.area(:,it)=Area.value(it,2*Lpf1-1:-2:1)';
+ gdat_data.dareadpsi(:,it)=Area.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
+ gdat_data.Rmesh(:,it) = Ri.value(it,1:M_Rmesh);
+ gdat_data.Zmesh(:,it) = Zj.value(it,1:N_Zmesh);
+ gdat_data.psi2D(1:M_Rmesh,1:N_Zmesh,it) = PFM_tree.value(1:M_Rmesh,1:N_Zmesh,it);
+ gdat_data.pressure(:,it)=Pres.value(it,2*Lpf1-1:-2:1)';
+ gdat_data.dpressuredpsi(:,it)=Pres.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
+ if ~isempty(Jpol.value)
+ gdat_data.jpol(:,it)=Jpol.value(it,2*Lpf1-1:-2:1)';
+ gdat_data.djpolpsi(:,it)=Jpol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
+ else
+ gdat_data.jpol = [];
+ gdat_data.djpolpsi = [];
+ end
+ gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)';
+ gdat_data.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1);
+ gdat_data.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1);
+ gdat_data.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1);
+ if ~isempty(Rinv.value)
+ gdat_data.rinv(:,it) = Rinv.value(it,Lpf1:-1:1)';
+ else
+ gdat_data.rinv = [];
+ end
+ if ~isempty(R2inv.value)
+ gdat_data.r2inv(:,it) = R2inv.value(it,Lpf1:-1:1)';
+ else
+ gdat_data.r2inv = [];
+ end
+ if ~isempty(Bave.value)
+ gdat_data.bave(:,it) = Bave.value(it,Lpf1:-1:1)';
+ else
+ gdat_data.bave = [];
+ end
+ if ~isempty(B2ave.value)
+ gdat_data.b2ave(:,it) = B2ave.value(it,Lpf1:-1:1)';
+ else
+ gdat_data.b2ave = [];
+ end
+ if ~isempty(FTRA.value)
+ gdat_data.ftra(:,it) = FTRA.value(it,Lpf1:-1:1)';
+ else
+ gdat_data.ftra = [];
+ end
+ %
+ end
+ gdat_data.x = gdat_data.rhopolnorm;
+ %
+ [equil_Rcoil,e]=rdaAUG_eff(shot,DIAG,'Rcl',exp_name_eff);
+ gdat_data.Rcoils=equil_Rcoil.value;
+ [equil_Zcoil,e]=rdaAUG_eff(shot,DIAG,'Zcl',exp_name_eff);
+ gdat_data.Zcoils=equil_Zcoil.value;
+ % get time values
+ [equil_time,e]=rdaAUG_eff(shot,DIAG,'time',exp_name_eff);
+ gdat_data.t = equil_time.value(1:NTIME);
+ %
+ [IpiPSI,e]=rdaAUG_eff(shot,DIAG,'IpiPSI',exp_name_eff);
+ gdat_data.Ip = IpiPSI.value(1:NTIME);
+ %
+ gdat_data.data = gdat_data.qvalue; % put q in data
+ gdat_data.units=[]; % not applicable
+ gdat_data.data_fullpath = [DIAG ' from expname: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
+ gdat_data.cocos = 17; % should check FPP
+ gdat_data.dim{1} = gdat_data.x;
+ gdat_data.dim{2} = gdat_data.t;
+ gdat_data.dimunits{1} = '';
+ gdat_data.dimunits{2} = 's';
- case 'nerho'
+ case {'ne','te'}
+ exp_name_eff = gdat_data.gdat_params.exp_name;
+ % ne or Te from Thomson data on raw z mesh vs (z,t)
+ nodenameeff = [upper(data_request_eff(1)) 'e_c'];
+ node_child_nameeff = [upper(data_request_eff(1)) 'e_core'];
+ [a,error_status]=rdaAUG_eff(shot,'VTA',nodenameeff,exp_name_eff);
+ if isempty(a.data) || isempty(a.t) || error_status>0
+ if nverbose>=3;
+ a
+ disp(['with data_request= ' data_request_eff])
+ end
+ return
+ end
+ gdat_data.(lower(node_child_nameeff)).data = a.data;
+ gdat_data.(lower(node_child_nameeff)).t = a.t;
+ gdat_data.t = a.t;
+ if any(strcmp(fieldnames(a),'units'))
+ gdat_data.units=a.units;
+ end
+ [alow,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'elow_c'],exp_name_eff);
+ [aup,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'eupp_c'],exp_name_eff);
+ gdat_data.(lower(node_child_nameeff)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff)).data,gdat_data.(lower(node_child_nameeff)).data-alow.value);
+ [a,error_status]=rdaAUG_eff(shot,'VTA','R_core',exp_name_eff);
+ gdat_data.(lower(node_child_nameeff)).r = repmat(a.data,size(gdat_data.(lower(node_child_nameeff)).data,1),1);
+ [a,error_status]=rdaAUG_eff(shot,'VTA','Z_core',exp_name_eff);
+ gdat_data.(lower(node_child_nameeff)).z = repmat(a.data',1,size(gdat_data.(lower(node_child_nameeff)).data,2));
+ gdat_data.dimunits=[{'Z [m]'} ; {'time [s]'}];
+ gdat_data.data_fullpath=[data_request_eff ' from VTA/' upper(data_request_eff(1)) 'e_c and ' upper(data_request_eff(1)) 'e_e'];
+ nb_core = size(gdat_data.(lower(node_child_nameeff)).z,1);
+ gdat_data.data = [];
+ gdat_data.data(1:nb_core,:) = gdat_data.(lower(node_child_nameeff)).data(1:nb_core,:);
+ gdat_data.dim=[{gdat_data.(lower(node_child_nameeff)).z};{gdat_data.(lower(node_child_nameeff)).t}];
+ gdat_data.error_bar(1:nb_core,:) = gdat_data.(lower(node_child_nameeff)).error_bar(1:nb_core,:);
+
+ % add edge part
+ nodenameeff_e = [upper(data_request_eff(1)) 'e_e'];
+ node_child_nameeff_e = [upper(data_request_eff(1)) 'e_edge'];
+ [a,error_status]=rdaAUG_eff(shot,'VTA',nodenameeff_e,exp_name_eff);
+ gdat_data.(lower(node_child_nameeff_e)).data = a.data;
+ ij_edge_notok = find(a.data>5e21);
+ gdat_data.(lower(node_child_nameeff_e)).data(ij_edge_notok) = NaN;
+ gdat_data.(lower(node_child_nameeff_e)).t = a.t;
+ if ~isempty(a.data)
+ [a,error_status]=rdaAUG_eff(shot,'VTA','R_edge',exp_name_eff);
+ gdat_data.(lower(node_child_nameeff_e)).r = repmat(a.data,size(gdat_data.(lower(node_child_nameeff_e)).data,1),1);
+ [a,error_status]=rdaAUG_eff(shot,'VTA','Z_edge',exp_name_eff);
+ gdat_data.(lower(node_child_nameeff_e)).z = repmat(a.data',1,size(gdat_data.(lower(node_child_nameeff_e)).data,2));
+ nb_edge = size(gdat_data.(lower(node_child_nameeff_e)).z,1);
+ iaaa=iroundos(gdat_data.(lower(node_child_nameeff_e)).t,gdat_data.(lower(node_child_nameeff)).t);
+ gdat_data.data(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).data(1:nb_edge,iaaa);
+ gdat_data.dim{1}(nb_core+1:nb_core+nb_edge,:)=gdat_data.(lower(node_child_nameeff_e)).z(1:nb_edge,iaaa);
+ [alow,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'elow_e'],exp_name_eff);
+ [aup,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'eupp_e'],exp_name_eff);
+ gdat_data.(lower(node_child_nameeff_e)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff_e)).data, ...
+ gdat_data.(lower(node_child_nameeff_e)).data-alow.value);
+ gdat_data.error_bar(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).error_bar(1:nb_edge,iaaa);
+ else
+ gdat_data.(lower(node_child_nameeff_e)).data = [];
+ gdat_data.(lower(node_child_nameeff_e)).t = [];
+ gdat_data.(lower(node_child_nameeff_e)).error_bar = [];
+ gdat_data.(lower(node_child_nameeff_e)).r = [];
+ gdat_data.(lower(node_child_nameeff_e)).z = [];
+ end
+ gdat_data.x=gdat_data.dim{1};
+ case {'ne_rho', 'te_rho'}
+ params_eff = gdat_data.gdat_params;
+ params_eff.data_request=data_request_eff(1:2);
+ % get raw data
+ [gdat_data,params_kin,error_status]=gdat_aug(shot,params_eff);
+ if error_status>0
+ if nverbose>=3; disp(['problems with ' params_eff.data_request]); end
+ return
+ end
+ % add rho coordinates
+ params_eff.data_request='equil';
+ [equil,params_equil,error_status]=gdat_aug(shot,params_eff);
+ if error_status>0
+ if nverbose>=3; disp(['problems with ' params_eff.data_request]); end
+ return
+ end
+ gdat_data.gdat_params.equil = params_equil.equil;
+
+ % core
+ node_child_nameeff = [upper(data_request_eff(1)) 'e_core'];
+ inb_chord_core=size(gdat_data.(lower(node_child_nameeff)).r,1);
+ inb_time_core=size(gdat_data.(lower(node_child_nameeff)).r,2);
+ psi_out_core = NaN*ones(inb_chord_core,inb_time_core);
+ rhopsinorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
+ rhotornorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
+ rhovolnorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
+ % edge
+ node_child_nameeff_e = [upper(data_request_eff(1)) 'e_edge'];
+ inb_chord_edge=size(gdat_data.(lower(node_child_nameeff_e)).r,1);
+ inb_time_edge=size(gdat_data.(lower(node_child_nameeff_e)).r,2);
+ psi_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
+ rhopsinorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
+ rhotornorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
+ rhovolnorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
+ % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
+ time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-1)];
+ for itequil=1:length(time_equil)-1
+ rr=equil.Rmesh(:,itequil);
+ zz=equil.Zmesh(:,itequil);
+ psirz_in = equil.psi2D(:,:,itequil);
+ it_core_inequil = find(gdat_data.(lower(node_child_nameeff)).t>=time_equil(itequil) & gdat_data.(lower(node_child_nameeff)).t<=time_equil(itequil+1));
+ if ~isempty(it_core_inequil)
+ rout_core=gdat_data.(lower(node_child_nameeff)).r(:,it_core_inequil);
+ zout_core=gdat_data.(lower(node_child_nameeff)).z(:,it_core_inequil);
+ psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core);
+ psi_out_core(:,it_core_inequil) = reshape(psi_at_routzout,inb_chord_core,length(it_core_inequil));
+ rhopsinorm_out_core(:,it_core_inequil) = sqrt((psi_out_core(:,it_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
+ for it_cx=1:length(it_core_inequil)
+ rhotornorm_out_core(:,it_core_inequil(it_cx)) = ...
+ interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
+ rhovolnorm_out_core(:,it_core_inequil(it_cx)) = ...
+ interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
+ end
+ end
+ % edge
+ it_edge_inequil = find(gdat_data.(lower(node_child_nameeff_e)).t>=time_equil(itequil) & gdat_data.(lower(node_child_nameeff_e)).t<=time_equil(itequil+1));
+ if ~isempty(it_edge_inequil)
+ rout_edge=gdat_data.(lower(node_child_nameeff_e)).r(:,it_edge_inequil);
+ zout_edge=gdat_data.(lower(node_child_nameeff_e)).z(:,it_edge_inequil);
+ psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_edge,zout_edge);
+ psi_out_edge(:,it_edge_inequil) = reshape(psi_at_routzout,inb_chord_edge,length(it_edge_inequil));
+ rhopsinorm_out_edge(:,it_edge_inequil) = sqrt((psi_out_edge(:,it_edge_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
+ for it_cx=1:length(it_edge_inequil)
+ rhotornorm_out_edge(:,it_edge_inequil(it_cx)) = ...
+ interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
+ rhovolnorm_out_edge(:,it_edge_inequil(it_cx)) = ...
+ interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
+ end
+ end
+ end
+ gdat_data.(lower(node_child_nameeff)).psi = psi_out_core;
+ gdat_data.(lower(node_child_nameeff)).rhopsinorm = rhopsinorm_out_core;
+ gdat_data.(lower(node_child_nameeff)).rhotornorm = rhotornorm_out_core;
+ gdat_data.(lower(node_child_nameeff)).rhovolnorm = rhovolnorm_out_core;
+ gdat_data.(lower(node_child_nameeff_e)).psi = psi_out_edge;
+ gdat_data.(lower(node_child_nameeff_e)).rhopsinorm = rhopsinorm_out_edge;
+ gdat_data.(lower(node_child_nameeff_e)).rhotornorm = rhotornorm_out_edge;
+ gdat_data.(lower(node_child_nameeff_e)).rhovolnorm = rhovolnorm_out_edge;
+ % put values of rhopsinorm for dim{1} by default
+ gdat_data.x = gdat_data.(lower(node_child_nameeff)).rhopsinorm;
+ iaaa=iroundos(gdat_data.(lower(node_child_nameeff_e)).t,gdat_data.(lower(node_child_nameeff)).t);
+ gdat_data.x(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhopsinorm(:,iaaa);
+ gdat_data.dim{1} = gdat_data.x;
+ gdat_data.dimunits{1} = 'rhopsinorm';
+ gdat_data.data_fullpath = [gdat_data.data_fullpath ' projected on equilibrium ' gdat_data.gdat_params.equil];
otherwise
warning(['switchcase= ' data_request_eff ' not known in gdat_aug'])
@@ -417,8 +816,6 @@ else
return
end
-if ishot==shot; mdsclose; end
-
gdat_data.mapping_for_aug = mapping_for_aug;
error_status=0;
diff --git a/crpptbx/AUG/geteqdskAUG.m b/crpptbx/AUG/geteqdskAUG.m
index d96a944d..fc6a8a3c 100644
--- a/crpptbx/AUG/geteqdskAUG.m
+++ b/crpptbx/AUG/geteqdskAUG.m
@@ -1,32 +1,38 @@
-function [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(shot,time,NR,NZ,savedir,deltaz,varargin);
+function [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(shot,time,zshift,varargin);
%
-% [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(shot,time,NR,NZ,savedir,deltaz,varargin);
+% [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(shot,time,zshift,varargin);
+%
+% if shot is a structure assumes obtained from gdat(shot,'equil',...);
%
% you can then do:
-% write_eqdsk('fname',eqdskAUG,17); % EQI is COCOS=17 by default, use [17,11] if you want an ITER version
+% write_eqdsk('fname',eqdskAUG,17); % EQI is COCOS=17 by default
% write_eqdsk('fname',eqdskAUG,[17 11]); % if you want an ITER version with COCOS=11
%
+if ~exist('shot') || isempty(shot);
+ disp('requires a shot or equil structure from gdat(shot,''equil'',...) in geteqdskAUG')
+ return
+end
+
if ~exist('time');
time_eff = 2.0;
else
time_eff = time;
end
-if ~exist('savedir')
- savedir_eff = '.';
+if ~exist('zshift')
+ zshift_eff = 0.; % no shift
else
- savedir_eff = savedir;
+ zshift_eff = zshift;
end
-if ~exist('deltaz')
- deltaz_eff = NaN; % no shift
+if isnumeric(shot)
+ equil=gdat(shot,'equil');
else
- deltaz_eff = deltaz;
+ equil = shot;
+ shot = equil.shot;
end
-
-equil=gdat(shot,'equil');
[zz it]=min(abs(equil.t-time_eff));
equil_all_t = equil;
@@ -36,13 +42,8 @@ eqdsk.cocos=17;
eqdsk.nr = size(equil.Rmesh,1);
eqdsk.nz = size(equil.Zmesh,1);
eqdsk.rmesh = equil.Rmesh(:,it);
-eqdsk.zmesh = equil.Zmesh(:,it);
-eqdsk.p = equil.pressure(:,it);
-eqdsk.pprime = equil.dpressuredpsi(:,it);
-eqdsk.FFprime = equil.ffprime(:,it);
-eqdsk.q = equil.qvalue(:,it);
-eqdsk.psimesh = equil.psi(:,it);
-eqdsk.rhopsi = equil.rhopolnorm(:,it);
+eqdsk.zmesh = equil.Zmesh(:,it) - zshift_eff;
+eqdsk.zshift = zshift_eff;
eqdsk.psi = equil.psi2D(:,:,it);
eqdsk.psirz = equil.psi2D(:,:,it);
eqdsk.rboxlen = equil.Rmesh(end,it) - equil.Rmesh(1,it) ;
@@ -51,8 +52,19 @@ eqdsk.zboxlen = equil.Zmesh(end,it) - equil.Zmesh(1,it) ;
eqdsk.zmid = 0.5*(equil.Zmesh(end,it) + equil.Zmesh(1,it)) ;
eqdsk.psiaxis = equil.psi_axis(it);
eqdsk.psiedge = equil.psi_lcfs(it);
+% need psimesh ot be equidistant and same nb of points as nr
+eqdsk.psimesh=linspace(0,1,eqdsk.nr)';
+eqdsk.rhopsi = sqrt(eqdsk.psimesh);
+% psimesh assumed from 0 on axis (so psi_edge is delta_psi) to have correct d/dpsi but to have sign(dpsi) easily
+eqdsk.psimesh = (eqdsk.psiedge-eqdsk.psiaxis) .* eqdsk.psimesh;
+nrho = size(equil.pressure,1);
+eqdsk.p = interpos(equil.rhopolnorm(:,it),equil.pressure(:,it),eqdsk.rhopsi,-0.03,[1 2],[0 equil.pressure(end,it)]);
+eqdsk.pprime = interpos(equil.rhopolnorm(:,it),equil.dpressuredpsi(:,it),eqdsk.rhopsi,-0.03,[1 2],[0 equil.dpressuredpsi(end,it)]);
+eqdsk.FFprime = interpos(equil.rhopolnorm(:,it),equil.ffprime(:,it),eqdsk.rhopsi,-0.03,[1 2],[0 equil.ffprime(end,it)]);
+eqdsk.q = interpos(equil.rhopolnorm(:,it),equil.qvalue(:,it),eqdsk.rhopsi,-0.03,[1 2],[0 equil.qvalue(end,it)]);
+
eqdsk.ip = equil.Ip(it);
-eqdsk.stitle=['AUGD equil_cliste t=' num2str(time_eff)];
+eqdsk.stitle=['#' num2str(shot) ' t=' num2str(time_eff) ' from ' equil.gdat_params.exp_name '/' equil.gdat_params.equil];
eqdsk.ind1=1;
psisign = sign(eqdsk.psimesh(end)-eqdsk.psimesh(1));
@@ -70,7 +82,7 @@ rmag=gdat(shot,'rmag');
[zz itrmag]=min(abs(rmag.t-time_eff));
eqdsk.raxis = rmag.data(itrmag);
zmag=gdat(shot,'zmag');
-eqdsk.zaxis = zmag.data(itrmag);
+eqdsk.zaxis = zmag.data(itrmag) - eqdsk.zshift;
% get plasma boundary
figure
@@ -108,4 +120,6 @@ eqdsk.rlim=AUG_innerouterwall(:,1);
eqdsk.zlim=AUG_innerouterwall(:,2);
plot(eqdsk.rlim,eqdsk.zlim,'k-')
+title(eqdsk.stitle)
+
eqdskAUG = eqdsk;
diff --git a/crpptbx/AUG/loadAUGdata.m b/crpptbx/AUG/loadAUGdata.m
index 8978a6eb..0d557599 100644
--- a/crpptbx/AUG/loadAUGdata.m
+++ b/crpptbx/AUG/loadAUGdata.m
@@ -1072,7 +1072,7 @@ switch AUGkeywrdcase{index}
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case {'te', 'ne'}
- shotfile_exp_eff = AUGexplocation{index};
+ shotfile_exp_eff = AUGexplocation{index}
if strcmp(AUGkeywrdcase{index},'te')
[a,e]=rdaAUG_eff(shot,'YPR','Te',shotfile_exp_eff);
diff --git a/crpptbx/AUG/rdaAUG_eff.m b/crpptbx/AUG/rdaAUG_eff.m
index 6feeeee3..6bf1e8ef 100644
--- a/crpptbx/AUG/rdaAUG_eff.m
+++ b/crpptbx/AUG/rdaAUG_eff.m
@@ -19,8 +19,7 @@ function [adata,error]=rdaAUG_eff(shot,diagname,sigtype,shotfile_exp,varargin);
%
global usemdsplus
-if isempty(usemdsplus); usemdsplus=1; end
-
+if ~exist('usemdsplus') || isempty(usemdsplus); usemdsplus=0; end
error=1;
time_int=[];
diff --git a/crpptbx/TCV/gdat_tcv.m b/crpptbx/TCV/gdat_tcv.m
index 56637e10..1093b7ad 100644
--- a/crpptbx/TCV/gdat_tcv.m
+++ b/crpptbx/TCV/gdat_tcv.m
@@ -32,7 +32,7 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(shot,data_req
% gdat_data.shot: shot number
% gdat_data.machine: machine providing the data
% gdat_data.gdat_request: keyword for gdat if relevant
-% gdat_data.data_fullpath: full path to the data node if known and relevant, or expression, or relevant function called if relevant
+% gdat_data.data_fullpath: full path to the data node if known and relevant, or relevant expression called if relevant
% gdat_data.gdat_params: copy gdat_params for completeness
% gdat_data.xxx: any other relevant information
%
@@ -69,7 +69,7 @@ gdat_params.doplot = 0;
gdat_params.liuqe = 1;
% construct list of keywords from global set of keywords and specific TCV set
-% get data_request names from centralized function
+% get data_request names from centralized table
data_request_names = get_data_request_names;
% add TCV specific to all:
if ~isempty(data_request_names.tcv)
@@ -354,24 +354,24 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_tcv.method,'expression')
% 2nd: method="expression"
- % assume expression contains function to call and which returns a structure into variable gdat_tmp
+ % assume expression contains an expression to evaluate and which creates a local structure into variable gdat_tmp
% we copy the structure, to make sure default nodes are defined and to avoid if return is an closed object like tdi
% eval_expr = [mapping_for_tcv.expression ';'];
eval([mapping_for_tcv.expression ';']);
if isempty(gdat_tmp) || (~isstruct(gdat_tmp) & ~isobject(gdat_tmp))
- warning(['function expression does not return a structure: ' mapping_for_tcv.expression])
+ warning(['expression does not create a gdat_tmp structure: ' mapping_for_tcv.expression])
error_status=801;
return
end
tmp_fieldnames = fieldnames(gdat_tmp);
if sum(strcmp(tmp_fieldnames,'data'))==0 % note: cannot do isfield since gdat_tmp might be an object
- warning(['function does not return a child name ''data'' for ' data_request_eff])
+ warning(['expression does not return a child name ''data'' for ' data_request_eff])
end
for i=1:length(tmp_fieldnames)
gdat_data.(tmp_fieldnames{i}) = gdat_tmp.(tmp_fieldnames{i});
end
% add .t and .x in case only dim is provided
- % do not allow shifting of timedim since should be treated in the relevant function
+ % do not allow shifting of timedim since should be treated in the relevant expression
ijdim=find(strcmp(tmp_fieldnames,'dim')==1);
if ~isempty(ijdim)
nbdims = length(gdat_data.dim);
@@ -397,7 +397,7 @@ elseif strcmp(mapping_for_tcv.method,'expression')
end
gdat_data.data_fullpath=mapping_for_tcv.expression;
end
- % end of method "function"
+ % end of method "expression"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_tcv.method,'switchcase')
@@ -592,8 +592,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
gdat_data.dim{1} = gdat_data.eqdsk{itime}.rmesh;
gdat_data.dim{2} = gdat_data.eqdsk{itime}.zmesh;
else
- aa=interpos2Dcartesian(gdat_data.eqdsk{itime}.rmesh,gdat_data.eqdsk{itime}.zmesh, ...
- gdat_data.eqdsk{itime}.psi,repmat(gdat_data.dim{1}',1,129),repmat(gdat_data.dim{2},129,1),-1,-1);
+ xx=repmat(reshape(gdat_data.dim{1},length(gdat_data.dim{1}),1),1,size(gdat_data.eqdsk{itime}.psi,2));
+ yy=repmat(reshape(gdat_data.dim{2},1,length(gdat_data.dim{2})),size(gdat_data.eqdsk{itime}.psi,1),1);
+ aa = interpos2Dcartesian(gdat_data.eqdsk{itime}.rmesh,gdat_data.eqdsk{itime}.zmesh ...
+ ,gdat_data.eqdsk{itime}.psi,xx,yy,-1,-1);
gdat_data.data(:,:,itime) = aa;
end
else
diff --git a/crpptbx/TCV/tcv_requests_mapping.m b/crpptbx/TCV/tcv_requests_mapping.m
index 94fe43f8..1cb4e89b 100644
--- a/crpptbx/TCV/tcv_requests_mapping.m
+++ b/crpptbx/TCV/tcv_requests_mapping.m
@@ -248,6 +248,7 @@ switch lower(data_request)
mapping.label = 'Zmagaxis';
mapping.method = 'tdiliuqe';
mapping.expression = '\results::z_axis';
+ %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% extra TCV cases (not necessarily in official data_request name list)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/crpptbx/gdat_plot.m b/crpptbx/gdat_plot.m
index 04a3c9b4..29974597 100644
--- a/crpptbx/gdat_plot.m
+++ b/crpptbx/gdat_plot.m
@@ -11,9 +11,8 @@ if ~isfield(gdat_data.gdat_params,'doplot') || gdat_data.gdat_params.doplot ==0
return
end
-fighandle = get(0,'CurrentFigure');
-
if prod(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(gdat_data.t)
+ fighandle = get(0,'CurrentFigure');
if gdat_data.gdat_params.doplot == 1
fighandle = figure;
elseif gdat_data.gdat_params.doplot > 1
@@ -25,7 +24,17 @@ if prod(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty
fighandle = figure(abs(gdat_data.gdat_params.doplot));
hold all
end
- if any(find(size(gdat_data.data)==length(gdat_data.t)))
+ if strcmp(gdat_data.gdat_request,'eqdsk')
+ % special case, plot contours of first equil
+ endstr = '';
+ if iscell(gdat_data.eqdsk); endstr = '{1}'; end
+ eval(['contour(gdat_data.eqdsk' endstr '.rmesh,gdat_data.eqdsk' endstr '.zmesh,gdat_data.eqdsk' endstr '.psi'',100);'])
+ hold on
+ eval(['plot(gdat_data.eqdsk' endstr '.rplas,gdat_data.eqdsk' endstr '.zplas,''k'');'])
+ eval(['plot(gdat_data.eqdsk' endstr '.rlim,gdat_data.eqdsk' endstr '.zlim,''k'');'])
+ axis equal;
+ title(eval(['gdat_data.eqdsk' endstr '.stitle']));
+ elseif 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')
@@ -33,7 +42,13 @@ if prod(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty
else
xlabel(['time']);
end
- ylabel([gdat_data.label '[' gdat_data.units ']']);
+ ylabel_eff = gdat_data.label;
+ if iscell(gdat_data.label) && length(gdat_data.label)>=2; ylabel_eff = gdat_data.label{2}; end
+ if ~isempty(gdat_data.units)
+ ylabel([ylabel_eff '[' gdat_data.units ']']);
+ else
+ ylabel(ylabel_eff);
+ end
zoom on;
end
else
--
GitLab