From 8fb3e50ecdd5dba2df1869a0a2cb442280083c71 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Thu, 17 Sep 2015 16:25:40 +0000
Subject: [PATCH] add and copy new gdat version
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@5023 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/AUG/aug_requests_mapping.m | 169 ++++
crpptbx/AUG/gdat_aug.m | 426 ++++++++
crpptbx/TCV/gdat_tcv.m | 1170 ++++++++++++++++++++++
crpptbx/TCV/tcv_requests_mapping.m | 275 +++++
crpptbx/gdat.m | 317 +++---
crpptbx/gdat_data_request_names_rho.xlsx | Bin 0 -> 12190 bytes
crpptbx/gdat_plot.m | 41 +
7 files changed, 2227 insertions(+), 171 deletions(-)
create mode 100644 crpptbx/AUG/aug_requests_mapping.m
create mode 100644 crpptbx/AUG/gdat_aug.m
create mode 100644 crpptbx/TCV/gdat_tcv.m
create mode 100644 crpptbx/TCV/tcv_requests_mapping.m
create mode 100644 crpptbx/gdat_data_request_names_rho.xlsx
create mode 100644 crpptbx/gdat_plot.m
diff --git a/crpptbx/AUG/aug_requests_mapping.m b/crpptbx/AUG/aug_requests_mapping.m
new file mode 100644
index 00000000..15b98dd6
--- /dev/null
+++ b/crpptbx/AUG/aug_requests_mapping.m
@@ -0,0 +1,169 @@
+function mapping = aug_requests_mapping(data_request)
+
+% 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)
+ 'min', -inf, ...
+ 'max', inf);
+
+if ~exist('data_request') || isempty(data_request)
+ return
+end
+
+% default label: data_request keyword itself
+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 = 'switchcase', then there will be a specific case within gdat_aug (usual case when not directly a signal)
+%
+% label is used for plotting
+switch lower(data_request)
+ case 'b0'
+ mapping.label = 'B_0';
+ mapping.method = 'switchcase';
+ mapping.expression = '';
+ case 'betan'
+ mapping.label = '\beta_N';
+ mapping.method = 'switchcase';
+ mapping.expression = '';
+ case 'betap'
+ mapping.label = '\beta_p';
+ mapping.method = 'tdi';
+ mapping.expression = '\results::beta_pol';
+ case 'cxrs'
+ 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'');'];
+ case 'delta_top'
+ mapping.label = 'delta\_top';
+ mapping.method = 'tdi';
+ mapping.expression = '\results::delta_ed_top';
+ case 'delta_bottom'
+ mapping.label = 'delta\_bottom';
+ mapping.method = 'tdi';
+ mapping.expression = '\results::delta_ed_bot';
+ case 'ece'
+ mapping.method = 'switchcase';
+ mapping.expression = '';
+ case 'eqdsk'
+ mapping.method = 'switchcase'; % could use function make_eqdsk directly?
+ mapping.expression = '';
+ case 'halpha'
+ mapping.label = 'Halpha';
+ mapping.method = 'tdi';
+ mapping.expression = '\base::pd:pd_011';
+ case 'ioh'
+ mapping.label = 'I ohmic transformer';
+ mapping.method = 'tdi';
+ mapping.expression = [{'\magnetics::ipol[*,$1]'} {'OH_001'}];
+ case 'ip'
+ mapping.label = 'Plasma current';
+ mapping.method = 'tdi';
+ mapping.expression = '\magnetics::iplasma:trapeze';
+ case 'kappa'
+ mapping.method = 'tdi';
+ mapping.expression = '\results::kappa_edge';
+ case 'mhd'
+ mapping.label = 'n=1,2, etc';
+ mapping.method = 'switchcase';
+ mapping.expression = '';
+ case 'ne'
+ mapping.method = 'switchcase';
+ case 'neint'
+ mapping.label = 'line integrated el. density';
+ mapping.method = 'tdi';
+ mapping.expression = '\results::fir:lin_int_dens';
+ case 'nel'
+ mapping.label = 'line-averaged el. density';
+ mapping.method = 'tdi';
+ mapping.expression = '\results::fir:n_average';
+ case 'nerho'
+ mapping.label = 'ne';
+ mapping.method = 'switchcase';
+ case 'neterho'
+ mapping.label = 'ne and Te';
+ mapping.method = 'switchcase';
+ case 'ni'
+ mapping.method = 'switchcase'; % especially since might have option fit, etc
+ case 'powers'
+ mapping.label = 'various powers';
+ mapping.method = 'switchcase';
+ case 'q0'
+ mapping.method = 'tdi';
+ mapping.expression = '\results::q_zero';
+ case 'q95'
+ mapping.method = 'tdi';
+ mapping.expression = '\results::q_95';
+ case 'qedge'
+ mapping.method = 'tdi';
+ mapping.expression = '\results::q_edge';
+ case 'qrho'
+ mapping.label = 'q';
+ mapping.method = 'switchcase';
+ case 'rgeom'
+ mapping.label = 'Rgeom';
+ mapping.method = 'switchcase';
+ 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';
+ case 'sxr'
+ mapping.method = 'switchcase';
+ case 'te'
+ mapping.label = 'Te';
+ mapping.method = 'switchcase';
+ case 'terho'
+ 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 = 'Volume';
+ mapping.method = 'switchcase';
+ % mapping.expression = '\results::psitbx:vol'; (if exists for liuqe2 and 3 as well)
+ case 'zeff'
+ mapping.label = 'zeff from Ip-Ibs';
+ mapping.method = 'tdi';
+ mapping.expression = '\results::ibs:z_eff';
+ case 'zgeom'
+ mapping.label = 'Zgeom';
+ mapping.method = 'switchcase';
+ case 'zmag'
+ mapping.label = 'Zmagaxis';
+ mapping.method = 'tdi';
+ mapping.expression = '\results::z_axis';
+% $$$ case ''
+% $$$ mapping.label = '';
+% $$$ mapping.method = 'tdi';
+% $$$ mapping.expression = '';
+ 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.expression = data_request;
+
+end
+
+
diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m
new file mode 100644
index 00000000..f070e344
--- /dev/null
+++ b/crpptbx/AUG/gdat_aug.m
@@ -0,0 +1,426 @@
+function [gdat_data,gdat_params,error_status,varargout] = gdat_aug(shot,data_request,varargin)
+%
+% function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
+%
+% Aim: get data from a given machine using full path or keywords.
+% data_request are and should be case independent (transformed in lower case in the function and outputs)
+%
+% If no inputs are provided, return the list of available pre-defined data_request in gdat_data and default parameters gdat_params
+%
+% Inputs:
+%
+% no inputs: return default parameters in a structure form in gdat_params
+% shot: shot number
+% data_request: keyword (like 'ip') or trace name or structure containing all parameters but shot number
+% varargin{i},varargin{i+1},i=1:nargin-2: additional optional parameters given in pairs: param_name, param_value
+% The optional parameters list might depend on the data_request
+% examples of optional parameters:
+% 'plot',1 (plot is set by default to 0)
+% 'machine','AUG' (the default machine is the local machine)
+%
+%
+% Outputs:
+%
+% gdat_data: structure containing the data, the time trace if known and other useful information
+% gdat_data.t : time trace
+% gdat_data.data: requested data values
+% gdat_data.dim : values of the various coordinates related to the dimensions of .data(:,:,...)
+% note that one of the dim is the time, replicated in .t for clarity
+% gdat_data.dimunits : units of the various dimensions, 'dimensionless' if dimensionless
+% gdat_data.error_bar : if provided
+% gdat_data.gdat_call : list of parameters provided in the gdat call (so can be reproduced)
+% 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.gdat_params: copy gdat_params for completeness
+% gdat_data.xxx: any other relevant information
+%
+%
+% Examples:
+% (should add working examples for various machines (provides working shot numbers for each machine...))
+%
+% [a1,a2]=gdat;
+% a2.data_request = 'Ip';
+% a3=gdat(48836,a2); % gives input parameters as a structure, allows to call the same for many shots
+% a4=gdat('opt1',123,'opt2',[1 2 3],'shot',48832,'data_request','Ip','opt3','aAdB'); % all in pairs
+% a5=gdat(48836,'ip'); % standard call
+% a6=gdat(48836,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (note all lowercase in output)
+
+%
+% Comments for local developer:
+% This gdat is just a "header routine" calling the gdat for the specific machine gdat_`machine`.m which can be called
+% directly, thus which should be able to treat the same type of input arguments
+%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Prepare some variables, etc
+
+varargout{1}=cell(1,1);
+error_status=1;
+
+% construct default parameters structure
+gdat_params.data_request = '';
+default_machine = 'aug';
+
+gdat_params.machine=default_machine;
+gdat_params.doplot = 0;
+
+% construct list of keywords from global set of keywords and specific AUG set
+% get data_request names from centralized function
+data_request_names = get_data_request_names;
+
+% add AUG specific to all:
+if ~isempty(data_request_names.aug)
+ aug_names = fieldnames(data_request_names.aug);
+ for i=1:length(aug_names)
+ data_request_names.all.(aug_names{i}) = data_request_names.aug.(aug_names{i});
+ end
+end
+data_request_names_all = fieldnames(data_request_names.all);
+
+% $$$ data_request_names_all= [{'ip'} {'b0'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'rhovol'} {'qrho'} {'q95'} {'kappa'} ...
+% $$$ {'delta'} {'deltatop'} {'deltabot'} {'neint'} {'nel'} ...
+% $$$ {'ne'} {'te'} {'nerho'} {'terho'} {'ne_edge'} {'te_edge'} {'nerho_edge'} {'terho_edge'} {'nerhozshift'} {'terhozshift'} {'profnerho'} {'profterho'} ...
+% $$$ {'neft'} {'teft'} {'neftav'} {'teftav'} {'neft:trial'} {'teft:trial'} {'neftav:trial'} {'teftav:trial'} ...
+% $$$ {'sxr'} {'sxr'} {'ece'} {'mpx'} {'ioh'} {'vloop'} {'pgyro'} {'jtor'} {'vi_tor'} {'vi_torfit'} {'vi_pol'} {'vi_polfit'} {'ti'} {'tifit'} {'ni'} {'nifit'} {'zeffcxrs'} {'zeffcxrsfit'}];
+
+% construct default output structure
+gdat_data.data = [];
+gdat_data.units = [];
+gdat_data.dim = [];
+gdat_data.dimunits = [];
+gdat_data.t = [];
+gdat_data.x = [];
+gdat_data.shot = [];
+gdat_data.gdat_request = [];
+gdat_data.gdat_params = gdat_params;
+gdat_data.data_fullpath = [];
+
+
+% Treat inputs:
+ivarargin_first_char = 3;
+data_request_eff = '';
+if nargin>=2 && ischar(data_request); data_request = lower(data_request); end
+
+gdat_data.gdat_request = data_request_names_all; % so if return early gives list of possible request names
+% no inputs
+if nargin==0
+ % return defaults and list of keywords
+ return
+end
+
+do_mdsopen_mdsclose = 1;
+% treat 1st arg
+if nargin>=1
+ if isempty(shot)
+ % means mdsopen(shot) already performed
+ shot = mdsipmex(2,'$SHOT');
+ gdat_data.shot = shot;
+ do_mdsopen_mdsclose = 0;
+ elseif isnumeric(shot)
+ gdat_data.shot = shot;
+ elseif ischar(shot)
+ ivarargin_first_char = 1;
+ else
+ warning('type of 1st argument unexpected, should be numeric or char')
+ error_status=2;
+ return
+ end
+ if nargin==1
+ % Only shot number given. If there is a default data_request set it and continue, otherwise return
+ return
+ end
+end
+% 2nd input argument if not part of pairs
+if nargin>=2 && ivarargin_first_char~=1
+ if isempty(data_request)
+ return
+ end
+ % 2nd arg can be a structure with all options except shot_number, or a string for the pathname or keyword, or the start of pairs string/value for the parameters
+ if isstruct(data_request)
+ if ~isfield(data_request,'data_request')
+ warning('expects field data_request in input parameters structure')
+ error_status=3;
+ return
+ end
+ data_request.data_request = lower(data_request.data_request);
+ data_request_eff = data_request.data_request;
+ gdat_params = data_request;
+ else
+ % since data_request is char check from nb of args if it is data_request or a start of pairs
+ if mod(nargin-1,2)==0
+ ivarargin_first_char = 2;
+ else
+ ivarargin_first_char = 3;
+ data_request_eff = data_request;
+ end
+ end
+end
+
+if ~isstruct(data_request)
+ gdat_params.data_request = data_request_eff;
+end
+
+% if start pairs from shot or data_request, shift varargin
+if ivarargin_first_char==1
+ varargin_eff{1} = shot;
+ varargin_eff{2} = data_request;
+ varargin_eff(3:nargin) = varargin(:);
+elseif ivarargin_first_char==2
+ varargin_eff{1} = data_request;
+ varargin_eff(2:nargin-1) = varargin(:);
+else
+ varargin_eff(1:nargin-2) = varargin(:);
+end
+
+% extract parameters from pairs of varargin:
+if (nargin>=ivarargin_first_char)
+ if mod(nargin-ivarargin_first_char+1,2)==0
+ for i=1:2:nargin-ivarargin_first_char+1
+ if ischar(varargin_eff{i})
+ % enforce lower case for any character driven input
+ if ischar(varargin_eff{i+1})
+ gdat_params.(lower(varargin_eff{i})) = lower(varargin_eff{i+1});
+ else
+ gdat_params.(lower(varargin_eff{i})) = varargin_eff{i+1};
+ end
+ else
+ warning(['input argument nb: ' num2str(i) ' is incorrect, expects a character string'])
+ error_status=401;
+ return
+ end
+ end
+ else
+ warning('number of input arguments incorrect, cannot make pairs of parameters')
+ error_status=402;
+ return
+ end
+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 ~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)
+ % copy description of keyword
+ gdat_data.request_description = data_request_names.all.(data_request_names_all{ij}).description;
+ end
+end
+
+% special treatment if shot and data_request given within pairs
+if isfield(gdat_params,'shot')
+ shot = gdat_params.shot; % should use only gdat_params.shot but change shot to make sure
+ gdat_data.shot = gdat_params.shot;
+ gdat_params=rmfield(gdat_params,'shot');
+end
+if ~isfield(gdat_params,'data_request') || isempty(gdat_params.data_request)
+ % warning('input for ''data_request'' missing from input arguments') % might be correct, asking for list of requests
+ error_status=5;
+ return
+end
+gdat_data.gdat_params = gdat_params;
+
+% re-assign main variables to make sure use the one in gdat_data structure
+shot = gdat_data.shot;
+data_request_eff = gdat_data.gdat_params.data_request;
+error_status = 6; % at least reached this level
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Specifications on how to get the data provided in aug_requests_mapping
+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
+ else
+ eval_expr = ['aatmp=tdi(''' mapping_for_aug.expression ''');'];
+ eval(eval_expr);
+ 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
+ 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)};
+ else
+ gdat_data.x = gdat_data.dim(dim_nontim);
+ end
+ end
+ gdat_data.t = gdat_data.dim{mapping_for_aug.timedim};
+ 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,...)
+ % 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)];
+ 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) ':)'];
+ % 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);
+ end
+ gdat_data.data_fullpath=mapping_for_aug.expression;
+ % end of method "tdi"
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+elseif strcmp(mapping_for_aug.method,'function')
+ % 2nd: method="function"
+ % assume expression contains function to call and which returns a structure
+ % 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])
+ 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])
+ end
+ for i=1:length(tmp_fieldnames)
+ gdat_data.(tmp_fieldnames{i}) = aatmp.(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
+ ijdim=find(strcmp(tmp_fieldnames,'dim')==1);
+ if ~isempty(ijdim)
+ 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
+ end
+ dim_nontim = setdiff([1:nbdims],mapping_for_aug.timedim);
+ ijt=find(strcmp(tmp_fieldnames,'t')==1);
+ if isempty(ijt)
+ gdat_data.t = gdat_data.dim{mapping_for_aug.timedim};
+ end
+ ijx=find(strcmp(tmp_fieldnames,'x')==1);
+ if isempty(ijx)
+ 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)};
+ else
+ gdat_data.x = gdat_data.dim(dim_nontim);
+ end
+ end
+ end
+ gdat_data.data_fullpath=mapping_for_aug.expression;
+ end
+ % end of method "function"
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+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])
+ return
+ end
+ if isempty(time) || ischar(time)
+ thomsontimes=time
+ warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check')
+ disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+ 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'];
+ 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;
+ end
+
+
+ case 'nerho'
+
+
+ otherwise
+ warning(['switchcase= ' data_request_eff ' not known in gdat_aug'])
+ error_status=901;
+ return
+ end
+
+else
+ warning(['AUG method=' mapping_for_aug.method ' not known yet, contact Olivier.Sauter@epfl.ch'])
+ error_status=602;
+ return
+end
+
+if ishot==shot; mdsclose; end
+
+gdat_data.mapping_for_aug = mapping_for_aug;
+error_status=0;
+
+return
+
diff --git a/crpptbx/TCV/gdat_tcv.m b/crpptbx/TCV/gdat_tcv.m
new file mode 100644
index 00000000..269a3fb2
--- /dev/null
+++ b/crpptbx/TCV/gdat_tcv.m
@@ -0,0 +1,1170 @@
+function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(shot,data_request,varargin)
+%
+% function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
+%
+% Aim: get data from a given machine using full path or keywords.
+% data_request are and should be case independent (transformed in lower case in the function and outputs)
+%
+% If no inputs are provided, return the list of available pre-defined data_request in gdat_data and default parameters gdat_params
+%
+% Inputs:
+%
+% no inputs: return default parameters in a structure form in gdat_params
+% shot: shot number
+% data_request: keyword (like 'ip') or trace name or structure containing all parameters but shot number
+% varargin{i},varargin{i+1},i=1:nargin-2: additional optional parameters given in pairs: param_name, param_value
+% The optional parameters list might depend on the data_request
+% examples of optional parameters:
+% 'plot',1 (plot is set by default to 0)
+% 'machine','TCV' (the default machine is the local machine)
+%
+%
+% Outputs:
+%
+% gdat_data: structure containing the data, the time trace if known and other useful information
+% gdat_data.t : time trace
+% gdat_data.data: requested data values
+% gdat_data.dim : values of the various coordinates related to the dimensions of .data(:,:,...)
+% note that one of the dim is the time, replicated in .t for clarity
+% gdat_data.dimunits : units of the various dimensions, 'dimensionless' if dimensionless
+% gdat_data.error_bar : if provided
+% gdat_data.gdat_call : list of parameters provided in the gdat call (so can be reproduced)
+% 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.gdat_params: copy gdat_params for completeness
+% gdat_data.xxx: any other relevant information
+%
+%
+% Examples:
+% (should add working examples for various machines (provides working shot numbers for each machine...))
+%
+% [a1,a2]=gdat;
+% a2.data_request = 'Ip';
+% a3=gdat(48836,a2); % gives input parameters as a structure, allows to call the same for many shots
+% a4=gdat('opt1',123,'opt2',[1 2 3],'shot',48832,'data_request','Ip','opt3','aAdB'); % all in pairs
+% a5=gdat(48836,'ip'); % standard call
+% a6=gdat(48836,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (note all lowercase in output)
+
+%
+% Comments for local developer:
+% This gdat is just a "header routine" calling the gdat for the specific machine gdat_`machine`.m which can be called
+% directly, thus which should be able to treat the same type of input arguments
+%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Prepare some variables, etc
+
+varargout{1}=cell(1,1);
+error_status=1;
+nverbose = 1;
+
+% construct default parameters structure
+gdat_params.data_request = '';
+default_machine = 'tcv';
+
+gdat_params.machine=default_machine;
+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
+data_request_names = get_data_request_names;
+% add TCV specific to all:
+if ~isempty(data_request_names.tcv)
+ tcv_names = fieldnames(data_request_names.tcv);
+ for i=1:length(tcv_names)
+ data_request_names.all.(tcv_names{i}) = data_request_names.tcv.(tcv_names{i});
+ end
+end
+data_request_names_all = fieldnames(data_request_names.all);
+
+% $$$ data_request_names_all= [{'ip'} {'b0'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'rhovol'} {'qrho'} {'q95'} {'kappa'} ...
+% $$$ {'delta'} {'deltatop'} {'deltabot'} {'neint'} {'nel'} ...
+% $$$ {'ne'} {'te'} {'nerho'} {'terho'} {'ne_edge'} {'te_edge'} {'nerho_edge'} {'terho_edge'} {'nerhozshift'} {'terhozshift'} {'profnerho'} {'profterho'} ...
+% $$$ {'neft'} {'teft'} {'neftav'} {'teftav'} {'neft:trial'} {'teft:trial'} {'neftav:trial'} {'teftav:trial'} ...
+% $$$ {'sxr'} {'sxr'} {'ece'} {'mpx'} {'ioh'} {'vloop'} {'pgyro'} {'jtor'} {'vi_tor'} {'vi_torfit'} {'vi_pol'} {'vi_polfit'} {'ti'} {'tifit'} {'ni'} {'nifit'} {'zeffcxrs'} {'zeffcxrsfit'}];
+
+% construct default output structure
+gdat_data.data = [];
+gdat_data.units = [];
+gdat_data.dim = [];
+gdat_data.dimunits = [];
+gdat_data.t = [];
+gdat_data.x = [];
+gdat_data.shot = [];
+gdat_data.gdat_request = [];
+gdat_data.gdat_params = gdat_params;
+gdat_data.data_fullpath = [];
+
+
+% Treat inputs:
+ivarargin_first_char = 3;
+data_request_eff = '';
+if nargin>=2 && ischar(data_request); data_request = lower(data_request); end
+
+gdat_data.gdat_request = data_request_names_all; % so if return early gives list of possible request names
+% no inputs
+if nargin==0
+ % return defaults and list of keywords
+ return
+end
+
+do_mdsopen_mdsclose = 1;
+% treat 1st arg
+if nargin>=1
+ if isempty(shot)
+ % means mdsopen(shot) already performed
+ shot = mdsipmex(2,'$SHOT');
+ gdat_data.shot = shot;
+ do_mdsopen_mdsclose = 0;
+ elseif isnumeric(shot)
+ gdat_data.shot = shot;
+ elseif ischar(shot)
+ ivarargin_first_char = 1;
+ else
+ warning('type of 1st argument unexpected, should be numeric or char')
+ error_status=2;
+ return
+ end
+ if nargin==1
+ % Only shot number given. If there is a default data_request set it and continue, otherwise return
+ return
+ end
+end
+% 2nd input argument if not part of pairs
+if nargin>=2 && ivarargin_first_char~=1
+ if isempty(data_request)
+ return
+ end
+ % 2nd arg can be a structure with all options except shot_number, or a string for the pathname or keyword, or the start of pairs string/value for the parameters
+ if isstruct(data_request)
+ if ~isfield(data_request,'data_request')
+ warning('expects field data_request in input parameters structure')
+ error_status=3;
+ return
+ end
+ data_request.data_request = lower(data_request.data_request);
+ data_request_eff = data_request.data_request;
+ gdat_params = data_request;
+ else
+ % since data_request is char check from nb of args if it is data_request or a start of pairs
+ if mod(nargin-1,2)==0
+ ivarargin_first_char = 2;
+ else
+ ivarargin_first_char = 3;
+ data_request_eff = data_request;
+ end
+ end
+end
+
+if ~isstruct(data_request)
+ gdat_params.data_request = data_request_eff;
+end
+
+% if start pairs from shot or data_request, shift varargin
+if ivarargin_first_char==1
+ varargin_eff{1} = shot;
+ varargin_eff{2} = data_request;
+ varargin_eff(3:nargin) = varargin(:);
+elseif ivarargin_first_char==2
+ varargin_eff{1} = data_request;
+ varargin_eff(2:nargin-1) = varargin(:);
+else
+ varargin_eff(1:nargin-2) = varargin(:);
+end
+
+% extract parameters from pairs of varargin:
+if (nargin>=ivarargin_first_char)
+ if mod(nargin-ivarargin_first_char+1,2)==0
+ for i=1:2:nargin-ivarargin_first_char+1
+ if ischar(varargin_eff{i})
+ % enforce lower case for any character driven input
+ if ischar(varargin_eff{i+1})
+ gdat_params.(lower(varargin_eff{i})) = lower(varargin_eff{i+1});
+ else
+ gdat_params.(lower(varargin_eff{i})) = varargin_eff{i+1};
+ end
+ else
+ warning(['input argument nb: ' num2str(i) ' is incorrect, expects a character string'])
+ error_status=401;
+ return
+ end
+ end
+ else
+ warning('number of input arguments incorrect, cannot make pairs of parameters')
+ error_status=402;
+ return
+ end
+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 ~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)
+ % copy description of keyword
+ gdat_data.request_description = data_request_names.all.(data_request_names_all{ij}).description;
+ end
+end
+
+% special treatment if shot and data_request given within pairs
+if isfield(gdat_params,'shot')
+ shot = gdat_params.shot; % should use only gdat_params.shot but change shot to make sure
+ gdat_data.shot = gdat_params.shot;
+ gdat_params=rmfield(gdat_params,'shot');
+end
+if ~isfield(gdat_params,'data_request') || isempty(gdat_params.data_request)
+ % warning('input for ''data_request'' missing from input arguments') % might be correct, asking for list of requests
+ error_status=5;
+ return
+end
+gdat_data.gdat_params = gdat_params;
+
+% re-assign main variables to make sure use the one in gdat_data structure
+shot = gdat_data.shot;
+data_request_eff = gdat_data.gdat_params.data_request;
+error_status = 6; % at least reached this level
+
+liuqe_version = 1;
+if isfield(gdat_data.gdat_params,'liuqe') && ~isempty(gdat_data.gdat_params.liuqe)
+ liuqe_version = gdat_data.gdat_params.liuqe;
+end
+substr_liuqe = '';
+if liuqe_version==2 || liuqe_version==3
+ substr_liuqe = ['_' num2str(liuqe_version)];
+end
+
+% special treatment for model shot=-1 or preparation shot >=100'000
+begstr = '';
+if shot==-1 || shot>=100000
+ % requires FBTE
+ liuqe_version = -1;
+ begstr = 'tcv_eq( "';
+ substr_liuqe = '", "FBTE" )';
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Specifications on how to get the data provided in tcv_requests_mapping
+mapping_for_tcv = tcv_requests_mapping(data_request_eff);
+gdat_data.label = mapping_for_tcv.label;
+
+ishot=NaN;
+if do_mdsopen_mdsclose
+ % mdsdefaultserver tcv1.epfl.ch; % should be in tcv general path, but set-it in the meantime...
+ if liuqe_version==-1
+ ishot = mdsopen('pcs', shot);
+ else
+ ishot = mdsopen(shot); % if ishot equal to shot, then mdsclose at the end
+ end
+ if ishot~=shot
+ warning(['cannot open shot= ' num2str(shot)])
+ return
+ end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% 1st treat the simplest method: "tdi" (and tdiliuqe)
+if strcmp(mapping_for_tcv.method(1:3),'tdi')
+ % need to treat liuqe2, model, etc from options....
+ substr_tdi = '';
+ if strcmp(mapping_for_tcv.method,'tdiliuqe'); substr_tdi = substr_liuqe; end
+ if iscell(mapping_for_tcv.expression)
+ if length(mapping_for_tcv.expression)>0
+ % series of arguments for tdi given in each cell
+ eval_expr = ['tdi(''' mapping_for_tcv.expression{1} substr_tdi ''''];
+ for i=2:length(mapping_for_tcv.expression)
+ eval_expr = [eval_expr ',''' mapping_for_tcv.expression{i} ''''];
+ end
+ eval_expr = [eval_expr ');'];
+ aatmp = eval(eval_expr);
+ else
+ % empty or wrong expression
+ error_status=701;
+ return
+ end
+ else
+ if liuqe_version==-1
+ mapping_for_tcv_expression_eff = mapping_for_tcv.expression;
+ if strcmp(lower(mapping_for_tcv.expression(1:8)),'\results')
+ mapping_for_tcv_expression_eff = mapping_for_tcv.expression(11:end);
+ end
+ eval_expr = ['tdi(''' begstr mapping_for_tcv_expression_eff substr_liuqe ''');']
+ else
+ eval_expr = ['tdi(''' mapping_for_tcv.expression substr_tdi ''');'];
+ end
+ aatmp=eval(eval_expr);
+ end
+ if isempty(aatmp.data) || isempty(aatmp.dim) % || ischar(aatmp.data) (to add?)
+ if (nverbose>=3); warning(['problems loading data for ' eval_expr ' for data_request= ' data_request_eff]); end
+ return
+ end
+ gdat_data.data = aatmp.data;
+ gdat_data.dim = aatmp.dim;
+ nbdims = length(gdat_data.dim);
+ if mapping_for_tcv.timedim==-1;
+ mapping_for_tcv.timedim = nbdims;
+ if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_tcv.timedim = nbdims-1; end
+ end
+ dim_nontim = setdiff([1:nbdims],mapping_for_tcv.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)};
+ else
+ gdat_data.x = gdat_data.dim(dim_nontim);
+ end
+ end
+ gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
+ gdat_data.units = aatmp.units;
+ gdat_data.dimunits = aatmp.dimunits;
+ 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.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.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 ';'];
+ eval(shift_eval);
+ 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];
+
+ % end of method "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
+ % 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])
+ 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])
+ 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
+ ijdim=find(strcmp(tmp_fieldnames,'dim')==1);
+ if ~isempty(ijdim)
+ nbdims = length(gdat_data.dim);
+ if mapping_for_tcv.timedim==-1;
+ mapping_for_tcv.timedim = nbdims;
+ if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_tcv.timedim = nbdims-1; end
+ end
+ dim_nontim = setdiff([1:nbdims],mapping_for_tcv.timedim);
+ ijt=find(strcmp(tmp_fieldnames,'t')==1);
+ if isempty(ijt)
+ gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
+ end
+ ijx=find(strcmp(tmp_fieldnames,'x')==1);
+ if isempty(ijx)
+ 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)};
+ else
+ gdat_data.x = gdat_data.dim(dim_nontim);
+ end
+ end
+ end
+ gdat_data.data_fullpath=mapping_for_tcv.expression;
+ end
+ % end of method "function"
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+elseif strcmp(mapping_for_tcv.method,'switchcase')
+ switch data_request_eff % not lower(...) since data_request_eff should be lower case already at this stage
+ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % First the request names valid for "all" machines:
+ %
+ case {'a_minor','rgeom'}
+ % compute average minor or major radius (on z=zaxis normally)
+ nodenameeff=['\results::r_max_psi' substr_liuqe];
+ rmaxpsi=tdi(nodenameeff);
+ nodenameeff2=['\results::r_min_psi' substr_liuqe];
+ rminpsi=tdi(nodenameeff2);
+ ij=find(rmaxpsi.data<0.5 | rmaxpsi.data>1.2);
+ if ~isempty(ij); rmaxpsi.data(ij)=NaN; end
+ ij=find(rminpsi.data<0.5 | rminpsi.data>1.2);
+ if ~isempty(ij); rminpsi.data(ij)=NaN; end
+ if strcmp(data_request_eff,'a_minor')
+ gdat_data.data=0.5.*(rmaxpsi.data(end,:) - rminpsi.data(end,:));
+ gdat_data.data_fullpath=[nodenameeff ' - ' nodenameeff2 ' /2'];
+ elseif strcmp(data_request_eff,'rgeom')
+ gdat_data.data=0.5.*(rmaxpsi.data(end,:) + rminpsi.data(end,:));
+ gdat_data.data_fullpath=[nodenameeff ' + ' nodenameeff2 ' /2'];
+ else
+ disp(['should not be in this case with data_request_eff = ' data_request_eff])
+ return
+ end
+ gdat_data.dim = rmaxpsi.dim(2);
+ gdat_data.t = gdat_data.dim{1};
+ if any(strcmp(fieldnames(rmaxpsi),'units'))
+ gdat_data.units = rmaxpsi.units;
+ end
+ gdat_data.dimunits = rmaxpsi.dimunits(2);
+
+ case {'zgeom'}
+ % compute average minor or major radius (on z=zaxis normally)
+ nodenameeff=['\results::z_contour' substr_liuqe];
+ zcontour=tdi(nodenameeff);
+ if strcmp(data_request_eff,'zgeom')
+ gdat_data.data=0.5.*(max(zcontour.data,[],1) + min(zcontour.data,[],1));
+ gdat_data.data_fullpath=['(max+min)/2 of ' nodenameeff];
+ gdat_data.dim{1} = zcontour.dim{2};
+ gdat_data.dimunits{1} = zcontour.dimunits{2};
+ else
+ disp(['should not be in this case with data_request_eff = ' data_request_eff])
+ return
+ end
+ gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
+ if any(strcmp(fieldnames(zcontour),'units'))
+ gdat_data.units = zcontour.units;
+ end
+
+ case {'b0'}
+ % B0 at R0=0.88
+ R0EXP=0.88;
+ if liuqe_version==-1
+ nodenameeff = 'tcv_eq("BZERO","FBTE")';
+ tracetdi=tdi(nodenameeff);
+ gdat_data.data = tracetdi.data;
+ else
+ nodenameeff=['\magnetics::iphi'];
+ tracetdi=tdi(nodenameeff);
+ gdat_data.data=192.E-07 * 0.996 *tracetdi.data/R0EXP;
+ end
+ if isempty(tracetdi.data) || isempty(tracetdi.dim) % || ischar(tracetdi.data) (to add?)
+ warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff])
+ return
+ end
+ gdat_data.data_fullpath=[nodenameeff];
+ gdat_data.dim = tracetdi.dim;
+ gdat_data.t = gdat_data.dim{1};
+ if any(strcmp(fieldnames(tracetdi),'units'))
+ gdat_data.units = tracetdi.units;
+ end
+ gdat_data.dimunits = tracetdi.dimunits;
+ gdat_data.request_description = ['vacuum magnetic field at R0=' num2str(R0EXP) 'm; COCOS=17'];
+
+ case {'betan'}
+ % 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;
+ params_eff.data_request='b0';
+ b0=gdat_tcv([],params_eff); % note: no need to set .doplot=0 since gdat_tcv does not call gdat_plot in any case
+ params_eff.data_request='ip';
+ ip=gdat_tcv([],params_eff);
+ params_eff.data_request='beta';
+ beta=gdat_tcv([],params_eff);
+ params_eff.data_request='a_minor';
+ a_minor=gdat_tcv([],params_eff);
+ % use beta as time base
+ if isempty(b0.data) || isempty(b0.dim) || isempty(ip.data) || isempty(ip.dim) || isempty(a_minor.data) || isempty(a_minor.dim) || isempty(beta.data) || isempty(beta.dim)
+ warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff])
+ return
+ 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{1} = 's';
+
+ 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;
+ else
+ disp(['"time" is expected as an option, choose default time = ' num2str(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;
+ for itime=1:length(time)
+ time_eff = time(itime);
+ % use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2
+ [fnames_readresults]=read_results_for_chease(shot,time_eff,liuqe_version,3,[],[],[],zshift,0,1);
+ eqdskval=read_eqdsk(fnames_readresults{4},7,0,[],[],1); % 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_eff,'%.4f')]);
+ % We still write COCOS=2 case, since closer to standard (in /tmp)
+ write_eqdsk(fnamefull,eqdsk_cocos_in,cocos_in,[],[],[],1);
+ % 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_eff (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]);
+ % 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) > 1
+ gdat_data.eqdsk{itime} = write_eqdsk(fnamefull,eqdsk_cocosout,cocos_out);
+ 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
+ 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);
+ gdat_data.data(:,:,itime) = aa;
+ end
+ else
+ 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;
+ 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 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 {'mhd'}
+ % load n=1, 2 and 3 Bdot from magnetic measurements
+ n1=tdi('abs(mhdmode("LFS",1,1))');
+ n2=tdi('abs(mhdmode("LFS",2,1))');
+ n3=tdi('abs(mhdmode("LFS",3,1))');
+ if ~isempty(n1.data)
+ gdat_data.data(:,1) = reshape(n1.data,length(n1.data),1);
+ if length(n2.data)==length(n1.data); gdat_data.data(:,2) = reshape(n2.data,length(n2.data),1); end
+ if length(n3.data)==length(n1.data); gdat_data.data(:,3) = reshape(n3.data,length(n3.data),1); end
+ gdat_data.dim{1} = n1.dim{1};
+ gdat_data.t = gdat_data.dim{1};
+ gdat_data.dim{2} = [1; 2; 3];
+ gdat_data.dimunits{1} = n1.dimunits{1};
+ gdat_data.dimunits{2} = 'n number';
+ gdat_data.units = 'T/s';
+ gdat_data.data_fullpath='abs(mhdmode("LFS",n,1))';
+ gdat_data.request_description = 'delta_Bdot from magnetic probes to get n=1, 2 and 3';
+ end
+
+ case {'ne','te'}
+ % ne or Te from Thomson data on raw z mesh vs (z,t)
+ edge_str_ = '';
+ edge_str_dot = '';
+ if isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
+ gdat_data.gdat_params.edge>0
+ edge_str_ = '_edge';
+ edge_str_dot = '.edge';
+ end
+ nodenameeff=['\results::thomson' edge_str_dot ':' data_request_eff];
+ tracetdi=tdi(nodenameeff);
+ tracestd=tdi(['\results::thomson' edge_str_dot ':' data_request_eff ':error_bar']);
+ 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])
+ return
+ end
+ if isempty(time) || ischar(time)
+ thomsontimes=time
+ warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check')
+ disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+ return
+ end
+ if strcmp(data_request_eff(1:2),'ne')
+ tracefirrat_data = get_fir_thom_rat_data(shot,['thomson' edge_str_],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'];
+ 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;
+ end
+
+ case {'ne_rho', 'te_rho', 'nete_rho'}
+ try
+ time=mdsdata('\results::thomson:times');
+ catch
+ warning('Problems with \results::thomson:times')
+ warning(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+ return
+ end
+ if isempty(time) || ischar(time)
+ thomsontimes=time
+ warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check')
+ disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+ return
+ end
+ 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
+ edge_str_ = '';
+ edge_str_dot = '';
+ if isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
+ gdat_data.gdat_params.edge>0
+ edge_str_ = '_edge';
+ edge_str_dot = '.edge';
+ end
+ % if nete_rho, do first ne, then Te later (so fir stuff already done)
+ if strcmp(data_request_eff,'ne_rho') || strcmp(data_request_eff,'nete_rho')
+ nodenameeff=['\results::thomson' edge_str_dot ':ne'];
+ tracetdi=tdi(nodenameeff);
+ nodenameeff=['\results::thomson' edge_str_dot ':ne; error_bar ; fir_thom_rat; (ne,std)*fir_thom_rat'];
+ tracestd=tdi(['\results::thomson' edge_str_dot ':ne:error_bar']);
+ tracefirrat_data = get_fir_thom_rat_data(shot,['thomson' edge_str_],time);
+ else
+ 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']);
+ end
+ gdat_data.data=tracetdi.data'; % Thomson data as (t,z)
+ gdat_data.error_bar=tracestd.data';
+ gdat_data.data_fullpath=nodenameeff;
+ if strcmp(data_request_eff,'ne_rho') || strcmp(data_request_eff,'nete_rho')
+ gdat_data.firrat=tracefirrat_data;
+ gdat_data.data_abs=gdat_data.data*diag(tracefirrat_data);
+ gdat_data.error_bar_abs=gdat_data.error_bar*diag(tracefirrat_data);
+ gdat_data.data_fullpath=[gdat_data.data_fullpath ' ; _abs includes *firrat'];
+ end
+ % add correct dimensions
+ % construct rho mesh
+ psi_max=tdi(['\results::thomson' edge_str_dot ':psi_max' substr_liuqe]);
+ psiscatvol=tdi(['\results::thomson' edge_str_dot ':psiscatvol' substr_liuqe]);
+ if abs(zshift)>1e-5
+ % calculate new psiscatvol
+ psitdi=tdi(['\results::psi' substr_liuqe]);
+ rmesh=psitdi.dim{1};
+ zmesh=psitdi.dim{2};
+ zthom=mdsdata('dim_of(\thomson:te,1)');
+ zeffshift=zshift;
+ % set zeffshift time array same as psitdi
+ switch length(zeffshift)
+ case 1
+ zeffshift=zeffshift * ones(size(psitdi.dim{3}));
+ case length(psitdi.dim{3})
+ % ok
+ case length(psiscatvol.dim{1})
+ zeffshift=interp1(psiscatvol.dim{1},zeffshift,psitdi.dim{3});
+ otherwise
+ disp(' bad time dimension for zshift')
+ disp(['it should be 1 or ' num2str(length(psiscatvol.dim{1})) ' or ' num2str(length(psitdi.dim{3}))])
+ end
+ for it=1:length(psiscatvol.dim{1})
+ itpsitdi=iround(psitdi.dim{3},psiscatvol.dim{1}(it));
+ psirz=psitdi.data(:,:,itpsitdi);
+ psiscatvol0=griddata(rmesh,zmesh,psirz',0.9*ones(size(zthom)),zthom-zeffshift(itpsitdi));
+ psiscatvol.data(it,:)=psiscatvol0;
+ end
+ end
+ if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data)
+ for ir=1:length(psiscatvol.dim{2})
+ rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))';
+ end
+ else
+ rho=NaN;
+ end
+ gdat_data.dim=[{rho};{time}];
+ gdat_data.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
+ gdat_data.x=rho;
+ gdat_data.t=time;
+ if any(strcmp(fieldnames(tracetdi),'units'))
+ gdat_data.units=tracetdi.units;
+ end
+ %%%%%%%%%%% add fitted profiles if 'fit'>=1
+ if isfield(gdat_data.gdat_params,'fit') && ~isempty(gdat_data.gdat_params.fit) && ...
+ gdat_data.gdat_params.fit>0
+ % default is from proffit:avg_time
+ def_proffit = '\results::proffit.avg_time';
+ if isfield(gdat_data.gdat_params,'fit_type') && ~isempty(gdat_data.gdat_params.fit_type)
+ if strcmp(gdat_data.gdat_params.fit_type,'local')
+ def_proffit = '\results::proffit.local_time';
+ else
+ gdat_data.gdat_params.fit_type = 'avg';
+ end
+ else
+ gdat_data.gdat_params.fit_type = 'avg';
+ end
+ if strcmp(data_request_eff(1:2),'ne')
+ nodenameeff = [def_proffit ':neft_abs']; % do first ne if nete asked for
+ elseif strcmp(data_request_eff(1:2),'te')
+ nodenameeff = [def_proffit ':teft'];
+ else
+ disp(['should not be here: data_request_eff, data_request_eff(1:2)= ',data_request_eff, data_request_eff(1:2)]);
+ end
+ if isfield(gdat_data.gdat_params,'trialindx') && ~isempty(gdat_data.gdat_params.trialindx) && ...
+ gdat_data.gdat_params.trialindx>=0
+ nodenameeff=[nodenameeff ':trial'];
+ trialindx = gdat_data.gdat_params.trialindx;
+ else
+ gdat_data.gdat_params.trialindx = [];
+ trialindx = [];
+ end
+ tracetdi=tdi(nodenameeff);
+ if isempty(trialindx)
+ gdat_data.fit.data = tracetdi.data;
+ else
+ if ~isempty(tracetdi.data) && size(tracetdi.data,3)>=trialindx+1
+ gdat_data.fit.data = tracetdi.data(:,:,trialindx+1);
+ else
+ gdat_data.fit.data = [];
+ gdat_data.fit.data_fullpath = [nodenameeff ' with trialindx=' num2str(trialindx) ' is empty'];
+ return
+ end
+ end
+ gdat_data.fit.x=tracetdi.dim{1};
+ gdat_data.fit.t=tracetdi.dim{2};
+ if mapping_for_tcv.timedim~=2 | mapping_for_tcv.gdat_timedim~=2
+ disp(['unexpected timedim in fit: data_request_eff= ' data_request_eff ...
+ ', mapping_for_tcv.timedim= ' mapping_for_tcv.timedim ...
+ ', mapping_for_tcv.gdat_timedim= ' mapping_for_tcv.gdat_timedim]);
+ end
+ gdat_data.dim=tracetdi.dim(1:2);
+ gdat_data.dimunits=tracetdi.dimunits(1:2);
+ if any(strcmp(fieldnames(tracetdi),'units'))
+ gdat_data.fit.units=tracetdi.units;
+ end
+ gdat_data.fit.data_fullpath = nodenameeff;
+ % do te as well if nete asked for
+ if strcmp(data_request_eff(1:4),'nete')
+ gdat_data.fit.ne.data = gdat_data.fit.data;
+ gdat_data.fit.ne.units = gdat_data.fit.units;
+ nodenameeff = [def_proffit ':teft'];
+ if ~isempty(trialindx); nodenameeff=[nodenameeff ':trial']; end
+ tracetdi=tdi(nodenameeff);
+ if isempty(trialindx)
+ gdat_data.fit.te.data = tracetdi.data;
+ else
+ if ~isempty(tracetdi.data) && size(tracetdi.data,3)>=trialindx+1
+ gdat_data.fit.te.data = tracetdi.data(:,:,trialindx+1);
+ else
+ return
+ end
+ end
+ if any(strcmp(fieldnames(tracetdi),'units'))
+ gdat_data.fit.te.units=tracetdi.units;
+ end
+ % construct pe=1.6022e-19*ne*te
+ gdat_data.fit.data = 1.6022e-19.*gdat_data.fit.ne.data .* gdat_data.fit.te.data;
+ gdat_data.fit.units = 'N/m^2; 1.6022e-19 ne Te';
+ gdat_data.fit.data_fullpath = [gdat_data.fit.data_fullpath ' ; ' nodenameeff ' and pe in data'];
+ end
+ else
+ gdat_data.gdat_params.fit = 0;
+ end
+ %%%%%%%%%%%
+ % if nete_rho, copy data as .ne, get .te and put pe=e ne Te in data:
+ if strcmp(data_request_eff(1:4),'nete')
+ gdat_data.ne.data = gdat_data.data_abs;
+ gdat_data.ne.error_bar = gdat_data.error_bar_abs;
+ gdat_data.ne.firrat=gdat_data.firrat;
+ gdat_data.ne.units = 'm^{-3}';
+ gdat_data = rmfield(gdat_data,{'firrat','data_abs','error_bar_abs'});
+ %
+ 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'];
+ 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 ...
+ + gdat_data.te.data .* gdat_data.ne.error_bar);
+ end
+
+ case {'powers'}
+ % note: same time array for all main, ec, ohm, nbi, ...
+ % At this stage fill just ech, later add nbi
+ nodenameeff='\results::toray.input:p_gyro';
+ tracetdi=tdi(nodenameeff);
+ gdat_data.ec.data = tracetdi.data*1e3; % at this stage p_gyro is in kW'
+ gdat_data.ec.units = 'W';
+ gdat_data.ec.dim=tracetdi.dim;
+ gdat_data.ec.dimunits=tracetdi.dimunits;
+ gdat_data.ec.t=tracetdi.dim{1};
+ gdat_data.ec.x=tracetdi.dim{2};
+ gdat_data.ec.data_fullpath=[nodenameeff];
+ gdat_data.ec.label='P_{EC}';
+ % set ec time as reference
+ gdat_data.t = gdat_data.ec.t;
+ gdat_data.dim{1} = gdat_data.t;
+ gdat_data.dimunits{1} = 's';
+ gdat_data.units = 'W';
+
+ % get ohmic power simply from vloop*Ip (minus sign for TCV)
+ ip=gdat([],'ip');
+ vloop=gdat([],'vloop');
+ tension = -1e5;
+ vloop_smooth=interpos(vloop.t,vloop.data,gdat_data.t,tension);
+ ip_t = interp1(ip.t,ip.data,gdat_data.t);
+ gdat_data.ohm.data = -vloop_smooth.*ip_t;
+ gdat_data.ohm.units = 'W';
+ gdat_data.ohm.dim=gdat_data.dim;
+ gdat_data.ohm.dimunits=gdat_data.dimunits;
+ gdat_data.ohm.t=gdat_data.t;
+ gdat_data.ohm.x=[];
+ gdat_data.ohm.data_fullpath=['-vloop(tens=' num2str(tension,'%.0e') ')*ip, from gdat'];
+ gdat_data.ohm.label='P_{OHM}';
+
+ % total power from each and total
+ gdat_data.data(:,1) = gdat_data.ohm.data;
+ gdat_data.data(:,2) = gdat_data.ec.data(:,10);
+ gdat_data.data(:,3) = gdat_data.ec.data(:,10) + gdat_data.ohm.data;
+ gdat_data.dim{2} = [1:3];
+ gdat_data.dimunits{2} = 'Pohm;Pec;Ptot';
+ gdat_data.data_fullpath=['tot power from EC and ohm'];
+ gdat_data.label = 'P_{ohm};P_{EC};P_{tot}';
+
+ case {'q_rho'}
+ % q profile on psi from liuqe
+ nodenameeff=['\results::q_psi' substr_liuqe];
+ if liuqe_version==-1
+ nodenameeff=[begstr 'q_psi' substr_liuqe];
+ end
+ tracetdi=tdi(nodenameeff);
+ gdat_data.data = tracetdi.data;
+ gdat_data.dim = tracetdi.dim;
+ gdat_data.t = gdat_data.dim{2};
+ gdat_data.data_fullpath=[nodenameeff ' on rhopol'];
+ rhopol_eff = ones(size(tracetdi.dim{1}));
+ rhopol_eff(:) = sqrt(linspace(0,1,length(tracetdi.dim{1})));
+ gdat_data.dim{1} = rhopol_eff;
+ gdat_data.x = gdat_data.dim{1};
+ gdat_data.dimunits{1} = '';
+ gdat_data.dimunits{2} = 's';
+ gdat_data.units = '';
+ gdat_data.request_description = 'q(rhopol\_norm)';
+
+ case {'psi_edge'}
+ % psi at edge, 0 by construction in Liuqe, thus not given
+ nodenameeff=['\results::psi_axis' substr_liuqe];
+ if liuqe_version==-1
+ nodenameeff=[begstr 'q_psi' substr_liuqe];
+ end
+ tracetdi=tdi(nodenameeff);
+ gdat_data.data = tracetdi.data.*0;
+ gdat_data.dim = tracetdi.dim;
+ gdat_data.t = gdat_data.dim{1};
+ gdat_data.data_fullpath=[' zero '];
+ gdat_data.dimunits = tracetdi.dimunits;
+ gdat_data.units = tracetdi.units;
+ gdat_data.request_description = '0 since LIUQE construct psi to be zero at LCFS';
+
+ case {'rhotor_edge','rhotor'}
+ % 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
+ % since cocos=17 for LIUQE we get:
+ % q = -dPhi/dpsi => Phi = - int(q*dpsi) which should always have the sign of B0
+ params_eff = gdat_data.gdat_params;
+ params_eff.data_request='q_rho';
+ q_rho=gdat_tcv([],params_eff);
+ params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE
+ psi_axis=gdat_tcv([],params_eff);
+ params_eff.data_request='b0'; % psi_edge=0 with LIUQE
+ b0=gdat_tcv([],params_eff);
+ b0tpsi = interp1(b0.t,b0.data,psi_axis.t); %q_rho on same time base as psi_axis
+ if isempty(psi_axis.data) || isempty(psi_axis.dim) || isempty(q_rho.data) || isempty(q_rho.dim)
+ warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff])
+ return
+ end
+ rhoequal = linspace(0,1,length(q_rho.dim{1}));
+ if strcmp(data_request,'rhotor_edge')
+ gdat_data.data = psi_axis.data; % to have the dimensions correct
+ gdat_data.dim = psi_axis.dim;
+ gdat_data.t = gdat_data.dim{1};
+ gdat_data.data_fullpath='phi from q_rho, psi_axis and integral(-q dpsi)';
+ gdat_data.units = 'T m^2';
+ gdat_data.dimunits{1} = 's';
+ elseif strcmp(data_request,'rhotor')
+ gdat_data.data = q_rho.data; % to have the dimensions correct
+ gdat_data.dim{1} = ones(size(q_rho.dim{1}));
+ gdat_data.dim{1}(:) = rhoequal;
+ gdat_data.dim{2} = q_rho.dim{2};
+ gdat_data.t = gdat_data.dim{2};
+ gdat_data.data_fullpath='sqrt(phitor/phitor_edge), rhotor_edge=sqrt(phitor/B0/pi)';
+ gdat_data.units = '';
+ gdat_data.dimunits{1} = 'rhopol\_norm';
+ gdat_data.dimunits{2} = 's';
+ end
+ for it=1:length(psi_axis.data)
+ ij=find(~isnan(q_rho.data(:,it)));
+ if ~isempty(ij)
+ [qfit,~,~,phi]=interpos(q_rho.x(ij).^2,q_rho.data(ij,it),rhoequal.^2);
+ dataeff = sqrt(phi .* psi_axis.data(it) ./ b0tpsi(it) ./ pi) ; % Delta_psi = -psi_axis
+ else
+ dataeff = NaN;
+ end
+ if strcmp(data_request,'rhotor_edge')
+ gdat_data.data(it) = dataeff(end);
+ elseif strcmp(data_request,'rhotor')
+ gdat_data.data(:,it) = dataeff./dataeff(end);
+ gdat_data.rhotor_edge(it) = dataeff(end);
+ end
+ gdat_data.b0 = b0tpsi(it);
+ end
+
+ case {'rhovol','volume_rho','volume'}
+ % volume_rho = vol(rho); volume = vol(LCFS) = vol(rho=1);
+ % rhovol = sqrt(vol(rho)/vol(rho=1));
+ nodenameeff='\results::psitbx:vol';
+ if liuqe_version==-1
+ nodenameeff=[begstr 'vol' substr_liuqe];
+ end
+ tracetdi=tdi(nodenameeff);
+ if isempty(tracetdi.data) || isempty(tracetdi.dim)
+ return
+ end
+ gdat_data.units = tracetdi.units;
+ if strcmp(data_request,'volume')
+ gdat_data.data = tracetdi.data(end,:);
+ gdat_data.dim{1} = tracetdi.dim{2};
+ gdat_data.data_fullpath=['\results::psitbx:vol(end,:)'];
+ gdat_data.dimunits{1} = tracetdi.dimunits{2};
+ gdat_data.request_description = 'volume(LCFS)=volume(rhopol=1)';
+ else
+ gdat_data.data = tracetdi.data;
+ gdat_data.dim = tracetdi.dim;
+ gdat_data.dimunits = tracetdi.dimunits;
+ if strcmp(data_request,'volume_rho')
+ gdat_data.data_fullpath=['\results::psitbx:vol'];
+ gdat_data.request_description = 'volume(rho)';
+ elseif strcmp(data_request,'rhovol')
+ gdat_data.volume_edge = gdat_data.data(end,:);
+ gdat_data.data = sqrt(gdat_data.data./repmat(reshape(gdat_data.volume_edge,1,size(gdat_data.data,2)),size(gdat_data.data,1),1));
+ gdat_data.data_fullpath='sqrt(\results::psitbx:vol/vol_edge)';
+ gdat_data.request_description = 'sqrt(volume(rho)/volume(edge))';
+ else
+ disp(['should not be here in vol cases with data_request = ' data_request_eff]);
+ return
+ end
+ end
+ gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
+
+ case {'sxr'}
+ % sxr from Xtomo by default or dmpx if 'camera','dmpx' is provided
+
+ case {'profnerho','profterho'}
+ % for backward compatibility but corresponds to ne_rho with param.fit_type='auto' (TCV special)
+ %
+ nodenameeff=['\results::THOMSON.PROFILES.AUTO:' data_request_eff(5:6)];
+ nodenameeff_vers = [nodenameeff ':version_num'];
+ avers = tdi(nodenameeff_vers);
+ if avers.data==0
+ % may be because nodes not yet filled in, so call once a node
+ ab=tdi(nodenameeff);
+ avers = tdi(nodenameeff_vers);
+ end
+ if avers.data>0
+ tracetdi=tdi(nodenameeff);
+ if avers.data < 2.99
+ % for earlier version the bug made it to have logically (rho,t)
+ gdat_data.data=tracetdi.data;
+ if ~isempty(tracetdi.dim) && ~ischar(tracetdi.data)
+ gdat_data.x=tracetdi.dim{1};
+ gdat_data.t=tracetdi.dim{2};
+ error_status=0;
+ else
+ error_status=2;
+ gdat_data.x=[];
+ gdat_data.t=[];
+ end
+ else
+ gdat_data.data=tracetdi.data'; % error in dimensions for autofits
+ if ~isempty(tracetdi.dim) && ~ischar(tracetdi.data)
+ if nverbose>=3; disp('assumes dim{2} for x in THOMSON.PROFILES.AUTO'); end
+ gdat_data.x=tracetdi.dim{2};
+ gdat_data.t=tracetdi.dim{1};
+ error_status=0;
+ else
+ gdat_data.x=[];
+ gdat_data.t=[];
+ error_status=2;
+ end
+ end
+ else
+ tracetdi=avers;
+ gdat_data.x=[];
+ gdat_data.t=[];
+ end
+ gdat_data.dim=[{gdat_data.x};{gdat_data.t}];
+ gdat_data.dimunits=[{'sqrt(psi\_norm)'} ; {'time [s]'}];
+ if ~isempty(gdat_data.t) && any(strcmp(fieldnames(tracetdi),'units'))
+ gdat_data.units=tracetdi.units;
+ end
+ gdat_data.request_description = 'quick autofits within thomson nodes, using version';
+ gdat_data.fullpath = ['Thomson autfits from ' nodenameeff];
+
+ otherwise
+ warning(['switchcase= ' data_request_eff ' not known in gdat_tcv'])
+ error_status=901;
+ return
+ end
+
+else
+ warning(['TCV method=' mapping_for_tcv.method ' not known yet, contact Olivier.Sauter@epfl.ch'])
+ error_status=602;
+ return
+end
+
+if ishot==shot; mdsclose; end
+
+gdat_data.mapping_for.tcv = mapping_for_tcv;
+error_status=0;
+
+return
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+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
+%
+% maintracename = 'thomson' or 'thomson_edge
+%
+% 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"
+%
+maintracename_eff = 'thomson';
+if exist('maintracename') && ~isempty(maintracename)
+ maintracename_eff = maintracename;
+end
+
+firthomratio = NaN;
+if ~exist('timebase') || isempty(timebase)
+ disp('need a timebase in get_fir_thom_rat_data')
+ 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)
+ disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat')
+ tracefirrat=tdi('\results::thomson:fir_thom_rat');
+ end
+ else
+ tracefirrat=tdi('\results::thomson:fir_thom_rat');
+ if isempty(tracefirrat.data) || ischar(tracefirrat.data)
+ disp('problem with \results::thomson:fir_thom_rat: empty')
+ 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)
+ disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat')
+ tracefirrat=tdi('\results::thomson:fir_thom_rat');
+ end
+ else
+ tracefirrat=tdi('\results::thomson_edge:fir_thom_rat');
+ if isempty(tracefirrat.data) || ischar(tracefirrat.data)
+ disp('problem with \results::thomson_edge:fir_thom_rat: empty')
+ end
+ end
+else
+ disp('bad input in get_fir_thom_rat_data')
+ return
+end
+
+if ~isempty(tracefirrat.data) || ischar(tracefirrat.data)
+ firthomratio = interp1(tracefirrat.dim{1},tracefirrat.data,timebase);
+end
+
diff --git a/crpptbx/TCV/tcv_requests_mapping.m b/crpptbx/TCV/tcv_requests_mapping.m
new file mode 100644
index 00000000..94fe43f8
--- /dev/null
+++ b/crpptbx/TCV/tcv_requests_mapping.m
@@ -0,0 +1,275 @@
+function mapping = tcv_requests_mapping(data_request)
+%
+% Information pre-defined for gdat_tcv, you can add cases here to match official cases in gdat_tcv, allowing backward compatibility
+%
+
+% Defaults
+mapping = struct(...
+ 'label', '', ...
+ 'method', '', ...
+ 'expression','', ...
+ '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
+end
+
+% default label: data_request keyword itself
+mapping.label = data_request;
+
+% for TCV, 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 = 'tdiliuqe': same as tdi but adds "_2" or "_3" if 'liuqe',2 or 3 is asked for in options
+% method = 'expression', then expression is executed and it should provide the structure gdat_tmp, which fields are copied to gdat_data
+% method = 'switchcase', then there will be a specific case within gdat_tcv (usual case when not directly a signal)
+%
+% 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 = 'expression';
+ % mapping.expression = ['tdi(''\results::q_psi'');']; % for tests
+ case 'delta_bottom'
+ mapping.timedim = 1;
+ mapping.label = 'delta\_bottom';
+ mapping.method = 'tdiliuqe';
+ mapping.expression = '\results::delta_ed_bot';
+ case 'delta_top'
+ mapping.timedim = 1;
+ mapping.label = 'delta\_top';
+ mapping.method = 'tdiliuqe';
+ mapping.expression = '\results::delta_ed_top';
+ 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 'ne_rho'
+ mapping.timedim = 2;
+ mapping.label = 'ne';
+ mapping.method = 'switchcase';
+ case 'neft'
+ mapping.timedim = 2;
+ mapping.label = 'ne';
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''ne_rho''; ' ...
+ 'params_eff.fit=1;params_eff.fit_type=''avg'';gdat_tmp=gdat_tcv([],params_eff);'];
+ case 'nete_rho'
+ 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.timedim = 1;
+ mapping.label = 'various powers';
+ mapping.method = 'switchcase';
+ case 'psi_axis'
+ mapping.timedim = 1;
+ mapping.method = 'tdiliuqe';
+ mapping.expression = '\results::psi_axis';
+ case 'psi_edge'
+ mapping.timedim = 1;
+ mapping.method = 'switchcase'; % is set to zero, so not in tree nodes
+ mapping.label = 'psi\_edge';
+ 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 'q_rho'
+ mapping.timedim = 2;
+ mapping.label = 'q';
+ mapping.method = 'switchcase';
+ case 'r_contour'
+ mapping.timedim = 1;
+ mapping.method = 'tdiliuqe';
+ mapping.expression = '\results::r_contour';
+ case 'rgeom'
+ mapping.timedim = 1;
+ mapping.label = 'Rgeom';
+ mapping.method = 'switchcase';
+ case 'rhotor_edge'
+ mapping.timedim = 1;
+ mapping.label = 'rhotor\_edge=sqrt(Phi/pi/B0)';
+ mapping.method = 'switchcase'; % from conf if exist otherwise computes it
+ case 'rhotor'
+ mapping.timedim = 2;
+ mapping.label = 'rhotor\_norm';
+ mapping.method = 'switchcase'; % from conf if exist otherwise computes it
+ 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 'te_rho'
+ mapping.timedim = 2;
+ mapping.label = 'Te';
+ mapping.method = 'switchcase';
+ case 'teft'
+ mapping.timedim = 2;
+ mapping.label = 'ne';
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''te_rho''; ' ...
+ 'params_eff.fit=1;params_eff.fit_type=''avg'';gdat_tmp=gdat_tcv([],params_eff);'];
+ 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 = 'Vloop';
+ mapping.method = 'tdi';
+ mapping.expression = [{'\magnetics::vloop[*,$1]'} {'001'}];
+ case 'volume'
+ mapping.timedim = 1;
+ mapping.label = 'Volume\_LCFS';
+ mapping.method = 'switchcase';
+ case 'volume_rho'
+ mapping.timedim = 2;
+ mapping.label = 'Volume(\rho)';
+ mapping.method = 'switchcase';
+ % mapping.expression = '\results::psitbx:vol'; (if exists for liuqe2 and 3 as well)
+ case 'z_contour'
+ mapping.timedim = 1;
+ mapping.method = 'tdiliuqe';
+ mapping.expression = '\results::z_contour';
+ 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';
+ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % extra TCV cases (not necessarily in official data_request name list)
+ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ case {'profnerho','profterho'}
+ mapping.timedim = 1;
+ mapping.label = data_request;
+ mapping.method = 'switchcase';
+ mapping.expression = '';
+
+% $$$ case ''
+% $$$ mapping.timedim = 1;
+% $$$ mapping.label = data_request;
+% $$$ mapping.method = 'tdi';
+% $$$ mapping.expression = '';
+ 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.expression = data_request;
+
+end
+
+if isempty(mapping.gdat_timedim)
+ mapping.gdat_timedim = mapping.timedim;
+end
diff --git a/crpptbx/gdat.m b/crpptbx/gdat.m
index 2b718b38..5b4804f9 100644
--- a/crpptbx/gdat.m
+++ b/crpptbx/gdat.m
@@ -1,206 +1,181 @@
-function [trace,error,varargout] = gdat(shot,data_type,varargin)
+function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
%
-% function [trace,error,varargout] = gdat(shot,data_type,varargin)
+% function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
%
-% new: no args or with shot=-9 to get keywords relevant for the related experiment
-% keywordlist=gdat; (for the default machine)
-% keywordlist=gdat(-9,[],'machine'); (for the specific machine)
+% Aim: get data from a given machine using full path or keywords.
+% Keywords should be the same for different machines, otherwise add "_machinname" to the keyword if specific
+% Keywords are and should be case independent (transformed in lower case in the function and outputs)
%
-% varargin{1}: 0 => no plot (default), 1 => plot (see special case for old JET style below)
-% varargin{2}: machine name 'AUG', 'JET' , 'TCV' (default: depending on where implemented)
+% If no inputs are provided, return the list of available keywords in gdat_data and default parameters gdat_params
%
-% list of data_type currently available:
+% Inputs:
%
-% (almost) All machines
-% 'Ip' = current
-% 'zmag' = vertical position of the center of the plasma (magnetic axis)
-% 'rmag' = radial position of the center of the plasma
-% 'rcont' = R of plama boundary vs time
-% 'zcont' = Z of plama boundary vs time
-% 'vol' = volume of flux surfaces vs rho=sqrt(psi)
-% 'qrho' = q profile on rho mesh
-% 'q95' = q95 vs time
-% 'kappa', 'elon' = edge elongation vs time
-% 'delta', 'triang' = edge averaged triangularity vs time
-% 'deltatop', 'triangtop' = edge upper (top) triangularity vs time
-% 'deltabot', 'triangbot' = edge lower (bottom) triangularity vs time
-% 'neint' = line-integrated electron density [m/m^3]
-% 'ne'= ne raw profile on (z,t). ADD error bars in .std
-% 'te'= Te raw profile on (z,t). ADD error bars in .std
-% 'nerho'= ne profile on (rho=sqrt(psi),time) mesh.Note rho is a 2D array as depends on time. ADD error bars in .std
-% 'terho'= Te profile on (rho=sqrt(psi),time) mesh.Note rho is a 2D array as depends on time. ADD error bars in .std
-% 'ece' = electron cyclotron emission
-% 'sxr' = soft x-ray emission
-% 'sxR' = soft x-ray emission with varargout{1} option (requires varargin{5}!)
-% 'Halpha' = H(D)-alpha trace
+% no inputs: return default parameters in a structure form in gdat_params
+% shot: shot number
+% data_request: keyword (like 'ip') or trace name or structure containing all parameters but shot number
+% varargin{i},varargin{i+1},i=1:nargin-2: additional optional parameters given in pairs: param_name, param_value
+% The optional parameters list might depend on the data_request
+% examples of optional parameters:
+% 'plot',1 (plot is set by default to 0)
+% 'machine','TCV' (the default machine is the local machine)
%
-% specific to TCV (see help loadTCVdata for more information)
-% 'xx_2 or xx_3' for Liuqe2 or 3: same as above for xx related to equilibrium
-% 'vol' = volume of flux surfaces
-% 'profnerho' = ne smoothed or fitted , vs (rho,t) (from Thomson fit)
-% 'profterho' = te smoothed or fitted , vs (rho,t) (from Thomson fit)
-% 'neft' = ne fitted from data on rho mesh (from proffit.local_time:neft)
-% 'teft' = te fitted from data on rho mesh (from proffit.local_time:teft)
-% 'neftav' = ne fitted from averaged over time data on rho mesh (from proffit.avg_time:neft)
-% 'teftav' = te fitted from averaged over time data on rho mesh (from proffit.avg_time:teft)
-% 'ece' = electron cyclotron emission
-% 'sxr' = soft x-ray emission
-% 'sxR' = soft x-ray emission with varargout{1} option (requires varargin{4}!)
-% 'MPX' = soft x-ray from wire chambers
-% 'vtor' = toroidal rotation, raw data and arror bars
-% 'vtorfit' = fitted toroidal rotation profiles (and error bars)
-% 'vpol' = poloidal rotation, raw data and arror bars
-% 'vpolfit' = fitted poloidal rotation profiles (and error bars)
-% 'Ti' (or Tc), 'Tifit', 'ni' (or nC), 'nifit', 'zeffcxrs', 'zeffcxrsfit': similar to 'vtor' from CXRS
%
-% JET
-% Special case compatible with old gdat.m allows (JET related): gdat(51994,'ppf','efit/xip',1)
+% Outputs:
%
-% INPUT:
+% gdat_data: structure containing the data, the time trace if known and other useful information
+% gdat_data.t : time trace
+% gdat_data.data: requested data values
+% gdat_data.dim : values of the various coordinates related to the dimensions of .data(:,:,...)
+% note that one of the dim is the time, replicated in .t for clarity
+% gdat_data.dimunits : units of the various dimensions, 'dimensionless' if dimensionless
+% gdat_data.error_bar : if provided
+% gdat_data.gdat_call : list of parameters provided in the gdat call (so can be reproduced)
+% 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.gdat_params: copy gdat_params for completeness
+% gdat_data.xxx: any other relevant information
%
-% shot: shot number
-% data_type: type of the required data.( see above)
-%
-% optional arguments valid for all values of data_type (not passed on to loadMACHINEdata function):
-%
-% Additional input arguments for specific traces (passed on to loadMACHINEdata function)
-%
-% data_type=sxr or ece:
-% varargin{3}: [i1 i2] : if not empty, assumes need many chords from i1 to i2
-% varargin{4}: channel status 1= unread yet, 0= read
-% (for traces with many channel, enables to load additional channels,
-% like SXR, ECE, etc.)
-% varargin{5}: zmag for varargout{1} computation
-%
-% OUTPUT:
-%
-% trace: structure containing data and meshes (compatible with tdi main structure)
-% trace.data: data
-% trace.t: time of reference
-% trace.x: space of reference
-% trace.dim: cell array of grids, trace.dim{1}, {2}, ...
-% trace.dimunits: units of dimensions
-% trace.std: if error bars are available
-%
-% error: error in loading signal (0=> OK, >=1 => error)
-%
-% Additional Output arguments depending on data_type
-%
-% data_type=sxR:
-% varargout{1}: intersection of the view lines with magnetic axis
-% data_type=MPX: (specific to TCV)
-% varargout{1}: see help loadTCVdata
-%
-%
-% functions needed: loadAUGdata, loadJETdata, loadTCVdata
-%
%
% Examples:
-% [zmag,error]=gdat(shot,'zmag',1); % gets zmag from TCV and plot
-% [zmag,error]=gdat(shot,'zmag',1,'JET'); % idem but from JET
-% [zmag,error]=gdat(shot,'ppf','efit/zmag',1); as above for JET
+% (should add working examples for various machines (provides working shot numbers for each machine...))
+%
+% [a1,a2]=gdat;
+% a2.data_request = 'Ip';
+% a3=gdat(48836,a2); % gives input parameters as a structure, allows to call the same for many shots
+% a4=gdat('opt1',123,'opt2',[1 2 3],'shot',48832,'data_request','Ip','opt3','aAdB'); % all in pairs
+% a5=gdat(48836,'ip'); % standard call
+% a6=gdat(48836,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (all lowercase in output)
+
+
%
-% keywordlist=gdat(-9,[],[],machine); % to get examples and the keywords defined for the relevant machine
-% keywordlist=gdat; % to get examples and the keywords defined for the local machine
+% Comments for local developer:
+% This gdat is just a "header routine" calling the gdat for the specific machine gdat_`machine`.m which can be called
+% directly, thus which should be able to treat the same type of input arguments
%
-gdatpaths
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Prepare some variables, etc
-nargineff=nargin;
-if nargineff>=3 & ischar(varargin{1})
-% $$$ if isempty(strmatch(upper(varargin{1}),[{'TCV'} {'JET'} {'AUG'}]))
- % old JET way with second part of data structure in varagin{1}
- data_type={data_type ; varargin{1}};
- if nargineff>=4;
- varargin{1}=varargin{2};
- else
- varargin{1}=0;
- nargineff=4;
+% for backward compatibility (with most recent simple call ala: gdat(shot,'data_request',doplot,machine) )
+% if nargin<=4 and >2 and 1st and 3rd are numeric, assume this is an old call and transform it to:
+% gdat(shot,'data_request','doplot',doplot','machine',machine)
+% and a warning
+if nargin>2
+ varargin_eff = varargin;
+ if nargin>=3 && nargin<=4 && isnumeric(shot) && isnumeric(varargin{1})
+ % assume old call: gdat(shot,'data_request',doplot[,machine])
+ varargin_eff{1} = 'doplot';
+ varargin_eff{end+1} = varargin{1};
+ if nargin==4
+ varargin_eff{end+1} = 'machine';
+ varargin_eff{end+1} = varargin{2};
+ end
end
- varargin{2}='JET';
-% $$$ end
-end
-
-% SETTING MACHINE
-doplot=0;
-if (nargineff>=3 & ~isempty(varargin{1})); doplot=varargin{1}; end
-a=which('gdat');
-if ~isempty(findstr('ipp',a)) | ~isempty(findstr('/u/osauter',a));
- machine='AUG';
- global usemdsplus; % so far from AUG, do not use mdsplus
- usemdsplus=0;
-elseif ~isempty(findstr('/home/osauter',a));
- machine='JET';
-elseif ~isempty(findstr('/home/sauter',a));
- machine='TCV';
-elseif ~isempty(findstr('/u/sauter',a));
- machine='D3D';
-end
-if ~exist('machine')
- % disp('did not find machine, set it to TCV')
- machine='TCV';
end
-if (nargineff>=4 & ~isempty(varargin{2})); machine=varargin{2}; end
-%machine
-% load data from specified machine
-if ~exist('data_type'); data_type=[]; end
-if ~exist('shot') || isempty(shot); shot=-9; end
-if nargineff<=4
- eval(['[trace,error,varargout] = load' machine 'data(shot,data_type);']);
+% construct default parameters structure
+gdat_params.data_request = '';
+fusion_machine_defaultname=getenv('FUSION_MACHINE_DEFAULTNAME');
+default_machine = '';
+if ~isempty(fusion_machine_defaultname)
+ default_machine = lower(fusion_machine_defaultname);
else
- eval(['[trace,error,varargout] = load' machine 'data(shot,data_type,varargin{3:end});']);
+ hostname=getenv('HOSTNAME');
+ if ~isempty(regexpi(hostname,'epfl'))
+ default_machine = 'tcv';
+ elseif ~isempty(regexpi(hostname,'rzg.mpg'))
+ default_machine = 'aug';
+ end
end
+machine_eff = default_machine;
-% PLOT DATA (if required)
-if doplot~=0; set_defaults_matlab; end
-if doplot==1 && isfield(trace,'data') && length(trace.data)>1 && ~ischar(trace.data)
- try
- figure;zoom on
- if length(size(trace.data))<=2
- hhh=plot(trace.t,trace.data);
- ylabel(data_type)
- else
- for idim=1:length(trace.dim)
- if length(trace.t)==length(trace.dim{idim}); idim_t=idim; end
- end
- if idim_t<=2
- hhh=plot(trace.t,trace.data(:,:,floor(end/2)));
- ylabel([data_type '(:,:,floor(end/2))'])
- elseif idim_t==3;
- hhh=plot(trace.t,reshape(trace.data(:,floor(end/2),:),length(trace.dim{1}),length(trace.t)));
- ylabel([data_type '(:,floor(end/2),:)'])
+% do not treat inputs here but in subsequent gdat_machine.m function, however, need to extract machine name if provided:
+if nargin>=2 % need at least 2 inputs to have 'machine','aug' as arguments (to ask specific list of requests)
+ % in case start directly with pairs of options
+ if ischar(shot) && ischar(data_request)
+ machine_eff = data_request;
+ elseif isstruct(shot) && isfield(shot,'machine')
+ machine_eff = shot.machine;
+ elseif isstruct(data_request) && isfield(data_request,'machine')
+ machine_eff = data_request.machine;
+ elseif nargin>=3
+ imachine=[];
+ for i=1:length(varargin_eff)
+ if ischar(varargin_eff{i})
+ ii = strmatch(varargin_eff{i},'machine','exact');
+ if ~isempty(ii); imachine = i; end
end
end
- xlabel('time [s]')
- title([machine ' ' num2str(shot)])
- grid on
- catch
- disp(' error in plotting part, most probably because could not guess time dimension correctly. To check')
+ if ~isempty(imachine) && length(varargin_eff)>=imachine+1 && ~isempty(varargin_eff{imachine+1})
+ machine_eff = varargin_eff{imachine+1};
+ end
+ end
+end
+
+% add paths to each machine which are in subdirectory of gdat and working
+% add only given machine directory using machine_eff...
+gdat_path = mfilename('fullpath');
+eval(['addpath ' gdat_path(1:end-4) upper(machine_eff)]);
+
+% copy gdat present call:
+if nargin==0
+ subcall=['gdat;'];
+elseif nargin>=1
+ if isnumeric(shot)
+ subcall=['gdat(' num2str(shot) ];
+ elseif ischar(shot)
+ subcall=['gdat(' shot ];
+ else
+ warning('type of 1st argument unexpected, should be numeric or char')
+ gdat_data.gdat_call = [];
+ return
end
-elseif doplot==-1
- try
- hold all
- child_h=get(gca,'child');
- nbplot=length(child_h);
- if length(size(trace.data))<=2
- hhh=plotos(trace.t,trace.data,'-',[],[],colos(mod(nbplot,12)+1,:));
+ if nargin>=2
+ if isempty(data_request)
+ subcall = [subcall ',[]'];
else
- hhh=plot(trace.t,trace.data(:,:,1),'--');
+ substring = subcall_all2str(data_request);
+ subcall = [subcall substring];
+ end
+ if nargin>=3
+ substring = subcall_all2str(varargin_eff{:});
+ subcall = [subcall substring];
end
- catch
- disp(' error in plotting part, most probably because could not guess time dimension correctly. To check')
end
+ % ... to continue
+ subcall = [subcall ');'];
end
+% Note: would need to check nargout to make call consistent, but to avoid this, each gdat_xxx should return at least an empty varargout: varargout{1}=cell(1);
try
- if exist('hhh') && ishandle(hhh(end)); set(hhh(end),'Tag',['gdat: ' num2str(shot)]); end
+ if nargin==0
+ eval(['[gdat_data,gdat_params,error_status,varargout] = gdat_' lower(machine_eff) ';']);
+ elseif nargin==1
+ eval(['[gdat_data,gdat_params,error_status,varargout] = gdat_' lower(machine_eff) '(shot);']);
+ elseif nargin==2
+ eval(['[gdat_data,gdat_params,error_status,varargout] = gdat_' lower(machine_eff) '(shot,data_request);']);
+ else
+ eval(['[gdat_data,gdat_params,error_status,varargout] = gdat_' lower(machine_eff) '(shot,data_request,varargin_eff{:});']);
+ end
catch
+ warning(['problems calling gdat_' lower(machine_eff)]);
+ return
end
-if doplot~=0
- h2=findobj(gca,'-regexp','Tag','gdat:*');
- if ~isempty(h2);
- legend(get(sort(h2),'Tag'));
+% copy subcall here so is last subnode
+gdat_data.gdat_call = [subcall ' % nargout = ' num2str(nargout)];
+
+if ~isfield(gdat_data.gdat_params,'doplot')
+ gdat_data.gdat_params.doplot = 0;
+end
+
+if gdat_data.gdat_params.doplot
+ % plot gdat_data versus 1st dim by default, if nb_dims<=2, otherwise do not plot
+ if length(varargout)==0 || isempty(varargout{1})
+ varargout{1} = gdat_plot(gdat_data); % return handle to figure
+ else
+ varargout{end+1} = gdat_plot(gdat_data); % return handle to figure
end
end
diff --git a/crpptbx/gdat_data_request_names_rho.xlsx b/crpptbx/gdat_data_request_names_rho.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..74f47e6263dadbe5c0879ab7e22d0c00441094c4
GIT binary patch
literal 12190
zcmeHtg<l+7()I**C%8+1V1v86y9Rf6nFNA60fI|#ch?|6g1fr~m*8&i<lf!)?q=`z
z{R8jz%&&WTPS?|Y&Z*~AbyX?JLPB8yU;uCc0Duf2%;`U7{|W$rf(8IE0dNqyq7L@1
zX7;WIY9L257d=K#J6qBmC<vM?00j8@|8D<>N1#k?z`lzatxfYlP^?`eEIXv61dgzU
zP?vHKHrX>)i?U%GkLJ<)jk&UVHC!i^vhA2I^#vbp`69Kph|mvP#(Ql=N<+lLPvy>M
z{7+Ju;{JJM8lP%oBG|KVkxx!wrs`qknih>2t;t`=sa~_!A+YXto-6E&l!gk5Qnk_T
zltA3gAMy`$X+y};y5Smgc?PI3;LEIPwzev>uxe0NrnwNQGsfS1qQ$CKYOcE{kv<rl
z#LCjEw#4HRB@wbJdsSq|LW0B^43Aa~#5yOn@P%QbSZNt-VnOI$^z<@3l*DaG|54H4
zX7I@CKM@nDo~1Zb<I!xr6^(e2;&O4s)}LP<#xSU3cBF6v|Llot=5^&jpLc=dH?~Kb
z0gO2VwFSbt1fyB409IR$(Ju*I$|3gP&n`7MaZYg_;!&cexwxXgECP}|jBq|z@LVB5
z`5_&NI$O?PM-twZE}^WX=LnOmC!V1YJcD`m@&XA^`Zt=^t20wvf^|<0+=>XMse!YZ
ztqT+5ul;}N`G1&`|I+ocSb4>6W~AT~smGB1o0*juG!a=3VX3cVY9D>2mr?6K<x&!?
zbi5%%Q^OB{l=Nx)_&B(<!uRQ*pX_>*tt<iqi<hFoqdYk2-pLJ?j@mg<!l`Vd2i<+<
zYUVmcQU*lh-WE+)+EkP+JFrF}F?A+ZjW))pi4TvNPZ<1~KgCa{S6+M7__pHJjEL%S
zd2m%DXVyX7SenmVV!=MLU?{iT(PRpCzl*WOLWNhqE!ou_zPhR<w^fyKhBFU2$iUS0
z+l5$K2j-I(gTj}76>>I=Tb41&KB}xMADue3v!QfQ4<~$gY5(bfe*{U<D%fKECrSJS
z;*&1GBhmsU2|fVsm8UJ!pFHt!aJDgaaIpDh&Hl-mS74I{M)|+Jm8r<f_JEgk#ODxZ
zk2H5|^aWQY@+0*_RJegU+9hf-){mEK1avKWTC=iDkPgA1laW668+MdUNX)Cww?*L?
zQ13BqPk5mnMlMDm;f#;ViTx#^F)@w~4oi;_QBqwnn+0O~=`*mnc5Z2?@IrD6$i*(3
zCxh9nN15GNtX;;52EjTrRJY}OAE<g%?qAo~NSaD(zWvfHx^xYb=?@~AxWI##=0gSn
zCibSu-rjbKU&EcZKAOveMr}HC?~%3%OJ#{m*!jJJT}eawB}sf}V_mj$QogF6Y`BU%
z4AOQiozo9po|tb9^v$e+XZk-?X4|Snbsek*{g40vD)<Vp%KWRc<i*KHPcftY;uP%A
zjh#S)UP93TLJR48;ZpT|eY*03(Y)$|XhYt-ryDs8gdqLbPEvdglF?7?KeU%g>W$YX
zKvwdjani$h_(rl;*W52H&#tVy!?4U$G&an;Z!RC(Js->`v_0mvKNvzW3+%eV$-frY
zpDyZDv<(P$r#obWN^SApA_!5s#|K(s=R0}fSj7ZK2M)Id8<ol+cva=hSa!iV6_Ji#
zc{33B6j+_%e}VC^Od~LI6AevQWBCF@W1(YO+E$I?9&Sn!w~!*YfvqFEQP!xQj7Epg
z1j^?O%caDl!ikWfTaIjCwM=EodQio(K=8WN8uT?AJ+ZXreKhH;9|_T*Ta+)DV%t$Z
zE5(YAO3)H6?Me-I>{Z=7;*g&bh!bnme3kC!D{BJ*x}AD@xw0{2vXPxJQoJosK-pKm
zY`{jm59l}CB_W$GbEDfuU2tz(w!NVWJPxXM)T{TYdM#E7B>((Fy!Y0j^Lvg@13_5&
zi=CdJS)YZK%+mIk)OA;MkM8qYdeQWYOOA}m=ca4C?RFII%nsSz?ggB^jqEA(?*3t=
zb^Q}W@6*kG3Oagr%{A20_Ag(W<y6)+%vipSLcNkvvU@O+`jTI3I?TmsbOxU8|5S!|
zuImx}U_TfU4FJIUQyE+=jhxL)Rb8E}>@8e=Eg$ikb`F@5XwUk^FGgFaG)f^Z<bDV}
z5|n`=&?@*yB}I!<_x3;S_uD}pV_bJv>20j*M~pEaPpMYhKt4MuX9H8^SXIJ@O&S*i
z6Hg0JO6n~;CFzQ*3<iYHJg;_BN>auVcH`c<#uY|n`tUf7LJrVzX1cDX^@(qkpisHo
z)C-GqerkvX-CHkp;YF0W6Xhi3aAk}Tzm9cc_mjC9j9Uqj#N5L+xbbY@<$lks{e`QY
z>n!3mwhpxrfuep(To;o<xFfrMrTJNyX}QSQTT8BR9uNs_#&UKTYG&I5gFs0Ghm$Wn
zdWS~tDZRiuol&(VMDym3a3;n28y)|bRo;Nv$-|zBd=APx`R#0_=yW4&>=?R@aeK#_
zZwO@95Oe&B{#KE+jtu)I`%~Xph%gHx0!GrdT4d$UDe)<i<B&>Lb<_knZkSki+-FDO
zYZwrdr6b!xcYL#@%mUmvS|Hg;cDf|ml#@DFodm7MuriU17RIJ7QJ%h8mT2~MR3sWX
z>uk@X?b7ssj5o>Fp0WJ7g37ghJd8Dkhk4bbMa0>dYc()W*lcCzi@k%h(ZoHHfmuCU
zT3fH{#UQ-Jb*LVS7^)Z)A8XTB_b73CnC$X6P`}D*xJ){*sROfwIZbUO(L3Kkpf<|q
zI4`MhG{|o(H&BY?T5$5+oj&4bIE6khg`U$OJEfl3okxvn+<t(VzGWOZk=P4-bH32r
zTOiGX=2-z3%>tQos8p|O_MCeDWv?Q|BA=>6IBC6BLQ`7+&HKwJGG?p75;Mv4`O(uZ
z2z_Y`5pC{8&~I%lX$?b*<`~_pkB=Yi`UmuPN*>4bpW?`Gp(NupNWCT-6bw5|gw~-J
z6ejr)2~WBV#d)D2&e8HQG}qbKO(GY(vLNAs_;+O^N=l0eHixK1ckjc#-&(AH{p#el
z`V7J&*<Vdt<+VIwo`;YTUDtP9B<g>NIiGO4r^dKO0Kl;94@_vW3nw<EIG~k8Exg`%
zy`ROu`F2j}v5xd8uB#JvprOMn=G_WIlmqn7o3JPn3`FELeJ=%R+;3uX<&#lf(jrK?
z$(8gcKQwe(`Md;%-w6%q%`Y2kvkRtn0k`ui{M;M;EwMR^gbwiiq$!$Cil<cK>K`~i
zNbzPJmkMdB7A5<ylHq21huL7HbiL8Ba(s-fB7L)kbaVEJ4c@H6r8fQ`ZG6O7T=!ks
zY1FCWaQdtu?V9K-*uq+$s0fpVv7cR?t217Kz$$q}>dtBeo6_~UX@bt4M@x>cIM1~D
zDd!`#AIjoy6fCw`IJ~oicIqfmMI2A>XTA!K8rcH0Zsb6ze0vx+@6^fFAow!<bH;+F
z#E^X=3xXQ6#ct~K^_3@&>B?C_-2)r*DJu}{%YOMAbG)~*&&aSzP-opk=$pb1hqX8k
zep(0bi9>lbSS!*>jezp~o`!*xh(DCps5~u5KZXT^%xxU~EEsw8$KQ-G=e~1ax+k~w
zbmEb%c?3m$g8X!8`fiX=-jAm)!#`nTV!u#2%_v3rb>%0+vS%q<@lf)~A-Hjegfj9$
zFW7K=jhU*R^i#H?TzT_vZ_bLkO@mkGM@2vb)Wb$nlSj|Fu(Td0l`aO*QJbD~&5jN_
z6ubBe7xa&Ll+$K8BYpQ&Rv#ipQ861kd~xO~etcr`hs+K0$8SYTm1^%lV*aX~=ybXX
z6H4n%aNk*Q(FxVE*jN@U0`uX+?1036r5sC7<PS6ZrF`+iLT64ProfW5dG$E`H&on{
zkqJj!WUVrnV%{Ag5IXu#n}UC9dwK0xWIMqvB4GPV@t5s&u{1Mtb@_+q_QwhEfBm+=
z=+5viW-Rd=$Y)W{4{@_WNTRBo<W1_kkjD3mFiS}>H&{;}Y6;E0ZErh{ZF8hO@D&r7
zqkN{HZh;8eGQqK|I{mQ7&>>vSb+C*D&%hUBG{xAS3aOwV6V?}_g^gGpq)hpQ+x@X;
zp*jk_I7VP+K{^b?Ja895FkV(-ZFg<}TgRkOe?q_Itr@_>ucG~Y6y;6AXP?_KDs?^<
zD9s&XD3Bw?-7kK8?rSW3=y9E4<BbAUU^hH;aSi`>28U)o?<fNX0Eou_6@2-_`*5{1
zvomA*WB*6DkF@n72sp94nJxv9E`HwHnA<X)!F;f_SURF6JIvcswizNeW2az^^)M01
zqKpm_vC_F$+0#YLw<Gmy&r^J&T`H)f-D0|wVZD~JX_lTs&ZDx?GG~+2l)TX=@=ot?
z<SiIQSKv01UMkNhZ9-9yu$7CoT51y@j8CoThL-6Kt$}J>f!_HUbnLOl9o--nXmTX<
z6zGuMKmZS1qD#7i3r831Gw1Siz)b9w`ewt!5UuKp<Zhx<p-v~vEkWs&4)xv89rx)G
z8TTtoN({PtW76ycR;m<K5;R5XPuV4|${eG<5cRNiMB>S5ZjriU5m6I&zMjUQL!};Q
zFUzmkSG$RgV@VSl`eEV7ZuRdWuaf45-bSL@l#IB(ud<%L@M794J9-PiK5Y7!OL7^I
zzz;2yi+;)WBWLOIYq|Z{ST#P)hr0Qn?C-Bi{N|44=D0k(?pRL|;umgaPXmPE^x43_
zH@qVh@$1Kb{k;*u?v*>a4?>2l6hfyO*uZbuguQZHSi&~#+Ka+oJ#Q1vO(bAaia{<t
zij|ksBInrFiL|eF1pwglq+)u|Bk(e^(DR%^X}G1rYGn)1DOZuGB1N6@HSp99P7R*<
zuJZW&8<;1pW%a)WR$@yFppfL7QZN2|Fr1ko9n18(c^shFy^9NZb3>fD{BXAQ{1o>=
z(6k{koR=o$06WwBX&}?!<KE-hUa1ADK<$s_sF$mU=#O{%s=R1Cx}F>RyQn@N9#8ji
zNM4?gArV{Sav+HWD1kepW?o0{On{g2pia{$)JY!=6I?RmX(9CahUA|S{Q!M^QtNKF
z>it#%TDVALPp6MV>n$YtF2eK-2I0n#O?|k>HPZ*avNa4N-Q6ITdd`edw`KOgC|kP_
z8lcEh=9rLn*oY~vQNWIRIEG&VS7+DfwRy7M#MKqbp3_5Oy!^=8Fh_fg!idlrDp6-c
zN)=_RkRG(8piiWMi>D*2$b~~i<#YUArsPIKq}F7DJ=2aq#My9Y8d3;zk`GYU#K~(u
zFB{}h5Y6pEkn?|!=Y9h2ENkq_AUw_Ze}wP(PJGSZ5-M*>*PlJI!)Pogc|C80o9$`7
zEf*`**WQOpCjQwa0iM7sI)y^4U6fF%2Dr=#$I<m<OxA?bQld(SA@(t)Fi}gtXP7V=
zkw;jo7^%NKdj}}3Pjg%-w&L^lV4=)aXPTF`8l^m4Riz#Nh6L-BIQ5MtFf!!|{;2$8
z$Be^`Pvu%#@A*QhJK1#B#&>JI@|7Q-WN+ZJT@C1IR=2IMMvdO%qh=}Dy`5flu1`o4
zKTUSPLRw}MFNum{IeUdZ@iEZ6s;9C<{xE*Gpl~?_i1OA)z3jT@Hr}q@C%Nx~?U!S%
zQuSuuk)8kw<#^3Ub?ze$v0xPZ?~Jx)`o}hf){U){q=(O#7&<1@VJ$&oCfabyx4al7
z78YGOW=0e7J4x9(5>DF)Sq<R`WtoIpZuNR=LkbXvszEfwLtJJd5JP)|>tC$$IJx5v
z@Bk5EW-wTv?WHRSzP|a8l)oq|wn&<tpC&CL+YR^0DI@)UVnDmMMLlpj>y2JDAu{wQ
zi3_(Lo$!8IyR5JZEH<a8C3~ft4Nu!Zb2|?RULVg$ix}Q_6(7*kyY3dG@)mr3u+38P
z-(8|bvxK~z7?>udhA2d#a#~7QF))Laj}NA?pTURgB9N_st~oFG9FUbJ+#3f^49|b0
z6L#r(jewA86?9cv@_b;Z9XCR@{4jAl&z*H~iSLhfL7=|>4Y=V5P#bCGrgqF}#;75G
zQIMBwK01&LrwgfZ7PCMfw7}9EW^KZlz$UtCIQP3wyul-E`Pmgv{KfF)s_bF2Sm(V0
z=7{kqo0|A&syoozJFd02Bv}kuBgM7U4w9*^qJEl3O6`#T%wPM4iW=&Z&P1}XJd1vk
zQSkz(So#<V2P*ZAKk`T25p)g3)`MZxD+P=4FQkcwG6}Z_KSevz4PhIF_})9B@D=)(
zDs`=hh8s>KrBhneyb|4N=${Kc-5OS%bU#shjKZ0HtWvYOe~Q(-q028=D^B8RH>kXR
zvd*6?UkM5m7jv2!`>6-xNGy(XcN{@y@YOFYNz#5-f}WAptzJc=X9iiW-hBC+=^X8-
z(rQ%`%BzWH%xReN_oW2!I#;7vTURY3DsPcm^7YX;-UsOKaq0$qt#HZvre8?>UbsUW
zp~Oh*+&?^T%{BpFzkK!#Kiz!oaCsyNIDR1~;B4Ug0lV_J<8KW*$pod=JMfzo6fyt+
z{{DKIc5wySn*H(O)T62Eu*i({F|+Ok?T&v0x(h&+B4d>&v%*%KUZ)wyWjK=54<vG3
z^jWKp6KZIPU*#YUW$IhI^BUd7eDmD2G2Bd6QlLUf!*%vqRE!akE)o9zeR{LW)--+u
z1S_qSe?y)G$<)osrj7$%vkrO=kuRE?zlM&c>R9Oo_C}MyXKZw4RhjSA!4*31-SXiW
zt!qPBnJ%9`_W`ksgp<!7anYh?ZrrwaVl}m<pv}}A7+yPj3bbh<6rcLPnvrMj0fpG{
zdaQGgB$SyID*{I2%%okMleT-(JbmMP^RZgoott0t2vhiS5ErZkJ(97oLcg*tpKNM(
zvu6%KCjnV~k)=+^O^kptOL6cq*D_SVuT>QXvWjLnO*Y*VY!M$~xC}NRi+NhyL7ckt
z9bLwe>Dx4a>bhtX1csspeTQib-c#vy4?~h4H6;BctdvUgS-8(?X>JK2k5cGC-nD#%
zVpcY&Q5g5?wdbr6K2sjudYKyrzVEI%%_lrs_B`2%?Lof!vUMM_dKjYuC*d;{gyBe`
z7VM`4M~MwCYN^U~E;r}g98v-mU$v0L2&aknT!X95)oM{Ju!BTDa*xtRo^IOsU+-|N
zW~F{%b~oF6EMLwaC%l#44`m-~YWLL=eZ~H!$;MxQ(LAoAz+`e>(x&R6>bcf2U4G$x
zxBB;n0-BW6t@%~jT-i=(Rw$$lg$K`Dw%&A33o7h3XzBv|o%IZw(q8#H<{&hGPlGMP
zw;}FawS9+vL>)00O&{H4U}H~T|J28Qm3uU&Fz5$|MC=H=$d&g>we<$UGm%^_-AgT}
z+lSy&qvMSTznBB&D_z`dtm}hgk{NeRP9m3dsc{RJPE&qm`q_9Zx}U{$zS_3e^-KD6
zb<OmPCG8gwB1}?Yb8-eJ)R=V+KIG7PF5^ZvS3#>%c4QcJ4hE^5AZ8X3asqGd?oynL
zQJo2>C%#PY<l(W}aiWLXYF@&a{&cb0W&-3f^9vBindr-r=E=`{|FqIJBhg~y;HlKk
zCyn|vy8(maK8)!@w##E6VA1<$oa0uegx`rdq@ZI61)qCguJD$*GTNwQ_W~$hod4&<
zMH6g04%Eo59~p@}`*xu6#CFM$k`Goa)Cy$UX^Wia!&uTK=ikz>vC!nzBiP^UV+R0m
z|D@rsnASi1>c2edGF@Gl%_fYGYb6&XG#yOqvKUf9wp388s%{E1(>w1nh?v7{0JROQ
zDZ`b548?YhEkQ<pNkXZLD+d~-SJDkXpKqJBZ$9RUr))2U?cRjwGp#V5n>Bc+d%oOx
zJU;XvCaBlwQ6GAxDYTB&T|PIx^f`EtfX>>yWT#4C@f)^#K9rWM92+Vzu1?CUTwmz5
z1Q}?6s-@{}8$F5Qzmm%*woI7w!v%LAc2W;CWK1~f66R7El-G#g$1|;R?l-W=guAAz
zjlY~TEG5)8jBMh%rRCi|w?BA)c5u(yq|&vRDbzT);)qzPtXyi9(w|+wQ52MUU2p#s
zt<*69j<W5PmY^QM6y#3so*D81<ZnENkoc60PO$f_H~DWH4RQz^9DNEQFW2awXa;Wd
z6I%HS4_%e>WpUqB7H+=8?yg~<z6}k|a!#_HS{>UZ^2`Qi9?#LOJ+bPkPr0^bct&gL
zGij?WL|+orC3`Y~)VI?FYv&qeqTQiu)IYz-G`0vXmK_vL3Z4mJi+n<vhWDm>nXKTu
zuQ2ZoJrJ{7X7hqe4L()B8#~f&%cUmqCzw?0vmb{#d9HZUhh>=^xGmpnqcr0gnR^*Z
zoFNcHi2ayS5qsk%9;hrB{%Gk*Bvf@eCwLLfR2+^`*L<+{4K?BhgM1E-Ojs>rnNcuv
z()krU!jDN<==B+L15wU(2Zy-g6d3r#M>Yzp_vlbPSrfW9s7q;yXpAcu0)lcwr!9i`
zdTY124U(gYWDS9ew)-JH=GxvRav>RRTG2|dUO7?UP9D0QJil@D$lTxy9O`$kyWzI<
z!klVzat)$c&IGU5x;ro9W%%LYrE@7tCoHdF-$|l|f0nWA_oT6VW>D+w;65N5?h8`6
z$hwr{!)5Fb#>f-=S()>-6IE&)TEwsmu0H`eg@<&GV@p~jn8SES$=O8I7}nHObO)t>
zKG+KwW0TMX31@++%08zUO*@>&QzSFd6Ch{r&tr`K!Bt4BrvL^bV9C;G!;=M$HDJba
zw25@bpBu*u2qCrAJLtyCtl{}Bb0|=fwx;_>ZG2p)uh|!2SLnjo?7=h?6r3K3@n8o9
zsKi|*MaCb^Dm@qhBmR=Gf#Md{q<ZZHP|`?@q0NTi_7f#crXOEbOdp7Al#@cyN!%g#
zomWUh58t7=bbf)OpEFW<*Xx$w@18!V6^2{ACx@Ij<Q70$;UYW??LuDlPIGyiJYR{n
zM$f|Rr!u(&J&!7+AJU30e10rYm4!8ig9nH+9QF+%c+u<~C9Vv5+OsW!6qIgj0n7w;
zgpQ=@Z4)U6gvxE`w{tHq+OY)k9d6xWVyYgCdZ%BHIVXd^)(pO2lNVj9{q3fAiqTx^
zynml2A5JiiV2|tFT9Fk5DJ4)w3a+RpsuRQ(r$oLM!#Eh`rT>BsIaQUu%Lk1%g$HS6
zJ-jzpwy!``IB(=NEHgt0<Z*WY!YD!H-Q8D)3n%jP8=S<)Wd0U8V6;UjnK~v__UbjQ
zK6@c9FCC6aWP1QhcoA+SHuvZkB~)<WvZ=6FjW>WxF;Z{E=>1Np)7Iek`8=lPoH0$!
z9xY{bw{7NZG4?=RIB9xH7rJWJAI4>qNqG~yZxM|xkt)Q^3C$Bf6?3pcHHsmXmp{zs
z`vBh;vsYgQ*+r1(;VLkoO^E{u*zJ(44;D34I$5d{Oc<3i;MZQBH`HCC$nXrmHEOg&
zQk<_d&)H7b0gJF3foNM8y4=P>``aL3%n{Jg<!OQwx0Lua9mbgF!dJB)k0a-8P*w_-
z$sD8%+gOBO-e<pcJQSog=i88GYLUeR$5|;I2+S%sG6|$|v<F;EA-p|p<%!fI1oYE{
zt{|fX*YM<uVasW@T2U$E_!wNJFmGx}ygl%}>X_HnS(tHc^4CkLj9YMNM5N5&q$4r#
zAx++Cg!<MDT^N%_6ZLu3PtP8`!>i?zN1wtASX4`coopIw^aNTL7?DXPx?3x~2kC23
zwh`8TKLax3CYO%v^pKpd!29#;!dk?!+KSN?#tawtpg@?#?LrIJ3<yC98;hU^blb8p
zqw$Q0G@y;4O^s&bQn9f|68A-BKc2T1aoy&Rah@E|5(-XSi?1+ejdGii-Wz3e*V?mj
z4nt6}Aw1xk3*qg6lw)?|+UMJCTpp7HDof!GBK?69xujgwCE3dJ9`VJQ-LrijyW-ge
zAUk+C?vBaXc_8PmX6Swy#uw$$D)e#;e-o8RJ*szvtVg0|Zh{#0#HP{`c<0_8>)%eN
z;Cfu&)8kTc<PalDNCO-C7$rcFeX#_={<L^V3<l3J4wU}gP#9$K)S*9Qdd9`pBymJ>
z^!U2WRz?*T2oxd-aeZhq28vPSS}U4Y_N+1oQzY_6DT15isAg$G7~MT1hoYwX$=0YX
zx63k+LHwdhg=6*AvlfhDL=7Sm!#km$OZE0PL;A#oF=`0fH*HqRGMqR<lPeZhimXgX
zGt?=9qaFkfa9WDifW*9^#FgthoIg5R{V58jOzs}F1$+|~vu^nPQS|suaHFYA2rgXE
z%hN!^LwLqK{1x)`Cc{*~l5CUK=<Xm!6c2E>B7?|FZF0m{-=mp$-HIuBC_`6VHF2z-
zviB%O-<;Ktf1iswn?Li?j;w~H=)FCqFkca0GU*&a6Xjp+8xavV;^O4(H5uq32Sy|J
zE)GQ;*fFQi=rRLWtQz#saS{lVL#gHofl=>2lItQ;VMcK2x)^ckipp82S*zOZ$o|+o
zLd2dLCXF1`Ur^j~$lh{zj_L5)hB}+8M4QHgR+uwa7$N6|3z%AJ<>rVhMoaJCcje`k
zq);%?!~M7(l3Y;(UPM$nvy&Z8OnPF}vM}Ca=U_@_*!LNK12Gtvk#`nC8e5W|h6KY|
zwWYxS$!UIp&*S2;<KkZXqLJG`MeiQhs1~l(ZVA(`)W0bIkLFK3xV+4d6F3L^ms)5o
zy`Hu6xq;?a0wOr+QhTnv?QhW;WoKb4g=;lwlvNAI{3qFWbTi@oSZ9!cQP%F{`4N=)
zUolSPd1%s$%tL7POiIelQGP^nf1eKrPhKn-<lL3_!*V6pws0hS<DEG&P;Urd#FG6-
zw5y3<ifR~tqar29aU(FZME{JR$Cy4^JcEEg9mFW5Z_+yA;<Ox$8@s!s5du=~7>+n-
zmMOui9EWlY>J8ECN$^T_NJ?qJo_@RrbodnuPs!*dgXCWluqqcCj><R8!&@p|dD9FN
zp*|7Wf3>8gOluBRdAWh{NKx_RwwMpOuN}N7d4{#6X&bY;Y7cnJ`nAq;^ke6Dge*$X
z=gxdIZ*8m1<E$j{SoKpm;;jznTcf8saQx=gyHDrS`DJvvS~;5>0cir?=klHMn*H5g
zeE!yB1RRn86~R$hBXHsz9h^8fbudwKc5rlIGI4M=`$q}bfAi(wp$m*PknLv1`o0F;
zAv*9OtzgL%(JsViK2EeVd{bBIXM<r}B4zDTBd1;CV)nYpPwpF!u{KIQVl_N3zv`Mv
zJaXkoi1Aho6$jQ!ow6=&gfbJIs+!4|YznQq;mS%CeN}xtPY!kDUg;{oH&(Ks^z3vO
z_V4rOL^S3UC*%6dY7)TRO!EkvUcrUeju?HMT7D+}+`+nrTBu^KA_xMchi$aHzU#O$
zN?GaODamj94!_YnScn!zevxT0IriZBPE-qh^1ai$%9)u%9+N$eF?dYz^x^(2Uza_;
zr41W4zR+LQ-nV&ht5HOnp`3NNb2ZPCvfhF|#m&a&>=Tre$8e^;%H!7qklrIP0{N+*
zCCN)!tNc%raO}aj<l!z);#d$mlV$9<=lKLm+nxjDHMAT46O$8qY#ym9Px_=M2YJRT
zP}0M>>ek=q;<=BYg-5`H90QJszxngM#Sxr$bv1KVHFI_Sm4IdX?HtHxY-R6kDrMAf
zA!RIOrep|b=xF3<LW;~>!(6puG~b2)*zWiFvtOLi>qdA<1V{=<!ETt>)TC{GL3mO;
zEXgsJi~`${O44W}P(927Eg8O~V^xv}YIV5Q@leUHm0*QoV9<hnk&zgNn`{H5ahc(n
zp^}+dEnTbD-ICLm9u>=kuu0Wl{Lqp}ich{OW5TM!W%WJ4Wc~2>7Ht)>$1>Dn-lG@-
zBn#`MzatcG#@fO$MRN9H7KIlj7bUL3k};)$h|M+J3I@>FfxG)R32FLoNw?8pUZR78
zDro=8OCv|e|I+f;nEknZiIsO)`c?V_c};?}LcqPkgXC9cya)SUs~RBvfk7)wA&@rq
z?Sw>)f#9Nqs>Mn%EF_k6hu5Xbfpf=&cnIxRTsmN8Le4kC95*A?u#l_c1rKFYf5Xl`
z-vlnh*GcymZ|dc$JH=33Ix)yJjb*3l^M=hMc>7ryVJoJA$G*4894xc-te!;jU%Q3Q
zTj@KbmoZ@?uFJR)!j+vSCM3n}*!T@Tz9YehnLUu4^0e^1O}t-C!Ez}pZdxF?lu?KX
z-`B8Q+}g{Mr7D!HXIVcB4{}2m$c=S?cMBO$Yk*3s(I8^rVZ&1Wyh74&{0?g@rn;|K
z40b@6wD(}Cw;Sz&T6MB<p66siGXffwje)}XhZbSizTnin_sm7MU%*taz9H0R7wy3!
zI_iTu6WN4M|E}t}pV1bv6UXblFhNh5+ezaBkD>rXzROIYpx$jdNvc~bOqV9FVw`n-
zP+%gbW8FitlE;zhd0YR?>n@mMo3RIxj>o>Upu@ZGWsfT7`fs?Y)>dnt5&x;guOJw}
z6%>D8WBkua`{(g*Dvp(8{~h4ps~`U<`0H2*X2!qNLjErJ?-gPHDmV;YWd7gvVZY=2
zUMBPx(hInZ`nRH?--UlKEBQ+p2|QrGf~3C}nf#9Md#3a+gb(0Sm4Dp+|CcxY9q{+W
z<6nT@=zjqIo`(D#;P+7PUjTC8kkhX|{cqs+chTR&Q-6t05&Yhnzk*f2qx`$e@)rgG
z08I@5{D%YcyZG<6<}ZK<nm-=mH%s%kmgjfSe_yWuf&>84Y5(P!{;_x~$wGs55CA{|
Of2hG~Y)t>_?0*2rA8rc(
literal 0
HcmV?d00001
diff --git a/crpptbx/gdat_plot.m b/crpptbx/gdat_plot.m
new file mode 100644
index 00000000..04a3c9b4
--- /dev/null
+++ b/crpptbx/gdat_plot.m
@@ -0,0 +1,41 @@
+function [fighandle]=gdat_plot(gdat_data,varargin);
+%
+% choices from doplot in gdat_data.gdat_params.doplot:
+% doplot = 0: no plot
+% = 1: new figure created
+% = -1: add to current figure (with hold all)
+% > 1: create new figure with this number, adding clf
+% <-1: add to figure number abs(doplot) (with hold all)
+%
+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)
+ if gdat_data.gdat_params.doplot == 1
+ fighandle = figure;
+ elseif gdat_data.gdat_params.doplot > 1
+ fighandle = figure(gdat_data.gdat_params.doplot);
+ clf;
+ elseif gdat_data.gdat_params.doplot == -1
+ hold all
+ elseif gdat_data.gdat_params.doplot < -1
+ fighandle = figure(abs(gdat_data.gdat_params.doplot));
+ hold all
+ end
+ 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
--
GitLab