diff --git a/crpptbx/JET/gdat_jet.m b/crpptbx/JET/gdat_jet.m
new file mode 100644
index 0000000000000000000000000000000000000000..5c3bcc7c25ff0467a3e9fcb1ec83272bb990843d
--- /dev/null
+++ b/crpptbx/JET/gdat_jet.m
@@ -0,0 +1,1507 @@
+function [gdat_data,gdat_params,error_status,varargout] = gdat_jet(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','JET' (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 relevant expression called if relevant
+% gdat_data.gdat_params: copy gdat_params for completeness (gdat_params contains a .help structure detailing each parameter)
+% gdat_data.xxx: any other relevant information
+%
+% gdat_params contains the options relevant for the called data_request. It also contains a help structure for each option
+% eg.: param1 in gdat_params.param1 and help in gdat_params.help.param1
+%
+% Examples:
+% (should add working examples for various machines (provides working shot numbers for each machine...))
+%
+% [a1,a2]=gdat;
+% a2.data_request = 'Ip';
+% a3=gdat(51994,a2); % gives input parameters as a structure, allows to call the same for many shots
+% a4=gdat('opt1',123,'opt2',[1 2 3],'shot',55379,'data_request','Ip','opt3','aAdB'); % all in pairs
+% a5=gdat(51994,'ip'); % standard call
+% a6=gdat(51994,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (note all lowercase in output)
+%
+% For JET, the specific trace can be given as:
+% aa==gdat(51994,[{'ppf'},{'efit'},{'xip'}],'machine','jet','doplot',1);
+%
+% 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 main default parameters structure
+gdat_params.data_request = '';
+default_machine = 'jet';
+
+gdat_params.machine=default_machine;
+gdat_params.doplot = 0;
+gdat_params.nverbose = 1;
+
+% construct list of keywords from global set of keywords and specific JET set
+% get data_request names from centralized table
+data_request_names = get_data_request_names;
+% add JET specific to all:
+if ~isempty(data_request_names.jet)
+ jet_names = fieldnames(data_request_names.jet);
+ for i=1:length(jet_names)
+ data_request_names.all.(jet_names{i}) = data_request_names.jet.(jet_names{i});
+ end
+end
+data_request_names_all = fieldnames(data_request_names.all);
+
+% 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 = [];
+gdat_data.help = [];
+
+% 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
+gdat_data.gdat_params.help = jet_help_parameters(fieldnames(gdat_data.gdat_params));
+gdat_params = gdat_data.gdat_params;
+
+% 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=-999; % means not defined
+ do_mdsopen_mdsclose = 0;
+ elseif isnumeric(shot)
+ gdat_data.shot = shot;
+ elseif ischar(shot)
+ ivarargin_first_char = 1;
+ else
+ if gdat_params.nverbose>=1; warning('type of 1st argument unexpected, should be numeric or char'); end
+ 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')
+ if gdat_params.nverbose>=1; warning('expects field data_request in input parameters structure'); end
+ 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
+ if gdat_params.nverbose>=1; warning(['input argument nb: ' num2str(i) ' is incorrect, expects a character string']); end
+ error_status=401;
+ return
+ end
+ end
+ else
+ if gdat_params.nverbose>=1; warning('number of input arguments incorrect, cannot make pairs of parameters'); end
+ 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:
+if ischar(data_request_eff) || length(data_request_eff)==1
+ ij=strmatch(data_request_eff,data_request_names_all,'exact');
+else
+ ij=[];
+end
+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 jet_requests_mapping
+mapping_for_jet = jet_requests_mapping(data_request_eff);
+gdat_data.label = mapping_for_jet.label;
+
+% fill again at end to have full case, but here to have present status in case of early return
+gdat_data.gdat_params.help = jet_help_parameters(fieldnames(gdat_data.gdat_params));
+gdat_data.mapping_for.jet = mapping_for_jet;
+gdat_params = gdat_data.gdat_params;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% 1st treat the simplest method: "signal"
+if strcmp(mapping_for_jet.method,'signal')
+ ppftype = mapping_for_jet.expression{1};
+ tracename = [mapping_for_jet.expression{2} '/' mapping_for_jet.expression{3}];
+ [a,x,t,d,e]=rda_jet(shot,ppftype,tracename);
+ if e==0
+ if gdat_params.nverbose>=3; disp(['error after rda_jet in signal with data_request_eff= ' data_request_eff]); end
+ return
+ end
+ aatmp.data = a; aatmp.t = t; aatmp.x = x;
+ gdat_data.data = aatmp.data;
+ gdat_data.t = aatmp.t;
+ gdat_data.x = aatmp.x;
+ if isempty(gdat_data.data)
+ return
+ end
+ if mapping_for_jet.timedim<=0
+ % need to guess timedim
+ if prod(size(aatmp.data)) == length(aatmp.data)
+ mapping_for_jet.timedim = 1;
+ elseif length(size(aatmp.data))==2
+ % 2 true dimensions
+ if length(aatmp.t) == size(aatmp.data,1)
+ mapping_for_jet.timedim = 1;
+ else
+ mapping_for_jet.timedim = 2;
+ end
+ else
+ % more than 2 dimensions, find the one with same length to define timedim
+ mapping_for_jet.timedim = find(size(aatmp.data)==length(aatmp.t));
+ if length(mapping_for_jet.timedim); mapping_for_jet.timedim = mapping_for_jet.timedim(1); end
+ end
+ mapping_for_jet.gdat_timedim = mapping_for_jet.timedim;
+ end
+ if length(size(aatmp.data))==1 | (length(size(aatmp.data))==2 & size(aatmp.data,2)==1)
+ gdat_data.dim=[{aatmp.t}];
+ gdat_data.dimunits={'time [s]'};
+ elseif length(size(aatmp.data))==2
+ gdat_data.dim{mapping_for_jet.timedim} = gdat_data.t;
+ gdat_data.dimunits{mapping_for_jet.timedim} = 'time [s]';
+ gdat_data.dim{setdiff([1:2],mapping_for_jet.timedim)} = gdat_data.x;
+ else
+ gdat_data.dim{mapping_for_jet.timedim} = gdat_data.t;
+ gdat_data.dimunits{mapping_for_jet.timedim} = 'time [s]';
+ nbdims = length(size(gdat_data.data));
+ for i=1:nbdims
+ if i~=mapping_for_jet.timedim
+ ieff=i;
+ if i>mapping_for_jet.timedim; ieff=i-1; end
+ gdat_data.dim{i} = gdat_data.x{ieff};
+ end
+ end
+ end
+ nbdims = length(gdat_data.dim);
+ if isfield(aatmp,'units')
+ gdat_data.units = aatmp.units;
+ elseif isfield(mapping_for_jet,'units')
+ gdat_data.units =mapping_for_jet.units;
+ end
+ keyboard
+ if mapping_for_jet.gdat_timedim>0 && mapping_for_jet.gdat_timedim ~= mapping_for_jet.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_jet.gdat_timedim-1);
+ inew=[1:mapping_for_jet.gdat_timedim-1 mapping_for_jet.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_jet.timedim-1) 'it,' ...
+ repmat(':,',1,nbdims-mapping_for_jet.timedim-1) ':)'];
+ dimstr_new=['(' repmat(':,',1,mapping_for_jet.gdat_timedim-1) 'it,' ...
+ repmat(':,',1,nbdims-mapping_for_jet.gdat_timedim-1) ':)'];
+ % eval gdat_data.data(;,:,...,it,...) = aatmp.data(:,:,:,it,...);
+ for it=1:size(aatmp.data,mapping_for_jet.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_jet.gdat_timedim = mapping_for_jet.timedim;
+ end
+ gdat_data.data_fullpath=mapping_for_jet.expression;
+ % end of method "signal"
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+elseif strcmp(mapping_for_jet.method,'expression')
+ % 2nd: method="expression"
+ % assume expression contains an expression to evaluate and which creates a local structure into variable gdat_tmp
+ % we copy the structure, to make sure default nodes are defined and to avoid if return is an closed object like tdi
+ % eval_expr = [mapping_for_jet.expression ';'];
+ eval([mapping_for_jet.expression ';']);
+ if isempty(gdat_tmp) || (~isstruct(gdat_tmp) & ~isobject(gdat_tmp))
+ if (gdat_params.nverbose>=1); warning(['expression does not create a gdat_tmp structure: ' mapping_for_jet.expression]); end
+ 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
+ if (gdat_params.nverbose>=1); warning(['expression does not return a child name ''data'' for ' data_request_eff]); end
+ 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 expression
+ ijdim=find(strcmp(tmp_fieldnames,'dim')==1);
+ if ~isempty(ijdim)
+ nbdims = length(gdat_data.dim);
+ if mapping_for_jet.timedim==-1;
+ mapping_for_jet.timedim = nbdims;
+ if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_jet.timedim = nbdims-1; end
+ end
+ dim_nontim = setdiff([1:nbdims],mapping_for_jet.timedim);
+ ijt=find(strcmp(tmp_fieldnames,'t')==1);
+ if isempty(ijt)
+ gdat_data.t = gdat_data.dim{mapping_for_jet.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_jet.expression;
+ end
+ % end of method "expression"
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+elseif strcmp(mapping_for_jet.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);
+ ijnan = find(isnan(rmaxpsi.data));
+ if isempty(rmaxpsi.data) || isempty(rmaxpsi.dim) || ischar(rmaxpsi.data) || ...
+ ( ~isempty(ijnan) && prod(size(ijnan))==prod(size(rmaxpsi.data)) )
+ if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
+ if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
+ return
+ end
+ nodenameeff2=['\results::r_min_psi' substr_liuqe];
+ rminpsi=tdi(nodenameeff2);
+ ijnan = find(isnan(rminpsi.data));
+ if isempty(rminpsi.data) || isempty(rminpsi.dim) || ...
+ ( ~isempty(ijnan) && prod(size(ijnan))==prod(size(rminpsi.data)) )
+ if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff2 ' for data_request= ' data_request_eff]); end
+ if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
+ return
+ end
+ 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
+ if (gdat_params.nverbose>=1); warning(['should not be in this case with data_request_eff = ' data_request_eff]); end
+ 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 isempty(zcontour.data) || isempty(zcontour.dim) % || ischar(zcontour.data) (to add?)
+ if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
+ if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
+ return
+ end
+ 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
+ if (gdat_params.nverbose>=1); warning(['should not be in this case with data_request_eff = ' data_request_eff]); end
+ return
+ end
+ gdat_data.t = gdat_data.dim{mapping_for_jet.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 = 'jet_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?)
+ if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
+ 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_jet, 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_jet([],params_eff); % note: no need to set .doplot=0 since gdat_jet does not call gdat_plot in any case
+ params_eff.data_request='ip';
+ ip=gdat_jet([],params_eff);
+ params_eff.data_request='beta';
+ beta=gdat_jet([],params_eff);
+ params_eff.data_request='a_minor';
+ a_minor=gdat_jet([],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)
+ if (gdat_params.nverbose>=1); warning(['problems loading data for data_request= ' data_request_eff]); end
+ 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_jet';
+ gdat_data.units = '';
+ gdat_data.dimunits{1} = 's';
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'cxrs'}
+ %not yet finished, just started
+ % 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 names from CXRS_get_profiles function, lower case when used in gdat
+ sub_nodes = {'Ti','vTor','vPol','nC','Zeff'}; % first node is also copied into data, choose "default' one
+ sub_nodes_out = lower({'Ti','vTor','vPol','ni','Zeff'}); %
+ sub_nodes_units = {'eV','m/s','m/s','m^{-3}',''}; % first node is also copied into data, choose "default' one
+ % use A. Karpushov routine to get profiles and then copy the data or the fitted profiles
+ aa=CXRS_get_profiles; cxrs_params = aa.param;
+ cxrs_params.k_plot=0; cxrs_params.k_debug=0;
+ % add params from gdat call
+ params_eff = gdat_data.gdat_params;
+ if isfield(params_eff,'cxrs_plot') && params_eff.cxrs_plot>0
+ cxrs_plot = params_eff.cxrs_plot;
+ else
+ cxrs_plot = 0;
+ end
+ gdat_data.gdat_params.cxrs_plot = cxrs_plot;
+ if isfield(params_eff,'source') && ~isempty(params_eff.source)
+ source = params_eff.source;
+ else
+ source = [1 2 3];
+ end
+ gdat_data.gdat_params.source = source;
+ if isfield(params_eff,'time_interval') && ~isempty(params_eff.time_interval) && length(params_eff.time_interval)>=2
+ time_interval = params_eff.time_interval;
+ cxrs_plot=1;
+ else
+ time_interval = [];
+ end
+ gdat_data.gdat_params.time_interval = time_interval;
+ gdat_data.gdat_params.cxrs_plot = cxrs_plot;
+ fit_tension_default = -100.;
+ if isfield(params_eff,'fit_tension')
+ fit_tension = params_eff.fit_tension;
+ else
+ fit_tension = fit_tension_default;
+ end
+ if ~isstruct(fit_tension)
+ fit_tension_eff.ti = fit_tension;
+ fit_tension_eff.vi = fit_tension;
+ fit_tension_eff.ni = fit_tension;
+ fit_tension_eff.zeff = fit_tension;
+ fit_tension = fit_tension_eff;
+ else
+ if ~isfield(fit_tension,'ti'); fit_tension.ti = fit_tension_default; end
+ if ~isfield(fit_tension,'vi'); fit_tension.vi = fit_tension_default; end
+ if ~isfield(fit_tension,'ni') && ~isfield(fit_tension,'nc'); fit_tension.ni = fit_tension_default; end
+ if ~isfield(fit_tension,'zeff'); fit_tension.zeff = fit_tension_default; end
+ end
+ gdat_data.gdat_params.fit_tension = fit_tension;
+ cxrs_params.prof.Ti.taus = fit_tension.ti;
+ cxrs_params.prof.vi.taus = fit_tension.vi;
+ cxrs_params.prof.nc.taus = fit_tension.ni;
+ cxrs_params.prof.zeff.taus = fit_tension.zeff;
+ cxrs_params.k_plot = cxrs_plot;
+ cxrs_profiles = CXRS_get_profiles(shot,source,time_interval,cxrs_params);
+ inb_times = length(cxrs_profiles.Times);
+ gdat_data.cxrs_params = cxrs_profiles.param;
+ if isempty(cxrs_profiles.Times) || ~isfield(cxrs_profiles,'proffit')
+ if (gdat_params.nverbose>=1); warning(['problems loading data with CXRS_get_profiles for data_request= ' data_request_eff]); end
+ for i=1:length(sub_nodes)
+ sub_eff_out = sub_nodes_out{i};
+ gdat_data.(sub_eff_out).fit.data = [];
+ gdat_data.(sub_eff_out).fit.rho = [];
+ gdat_data.(sub_eff_out).fit.error_bar = [];
+ gdat_data.(sub_eff_out).raw.data = [];
+ gdat_data.(sub_eff_out).raw.rho = [];
+ gdat_data.(sub_eff_out).raw.error_bar = [];
+ gdat_data.(sub_eff_out).raw.error_bar_rho = [];
+ gdat_data.(sub_eff_out).raw.cxrs_system = [];
+ gdat_data.time_interval = [];
+ gdat_data.(sub_eff_out).units = sub_nodes_units{i};
+ if i==1
+ gdat_data.data = gdat_data.(sub_eff_out).fit.data;
+ gdat_data.x = gdat_data.(sub_eff_out).fit.rho;
+ gdat_data.error_bar = gdat_data.(sub_eff_out).fit.error_bar;
+ gdat_data.units = gdat_data.(sub_eff_out).units;
+ end
+ end
+ return
+ end
+ inb_channels =120; % need to change if gets bigger!!! but easier to prefill with NaNs and use the "use" part
+ for i=1:length(sub_nodes)
+ sub_eff = sub_nodes{i};
+ sub_eff_out = sub_nodes_out{i};
+ % fits
+ if isfield(cxrs_profiles.proffit,sub_eff)
+ gdat_data.(sub_eff_out).fit.data = cxrs_profiles.proffit.(sub_eff);
+ gdat_data.(sub_eff_out).fit.rho = cxrs_profiles.proffit.([sub_eff '_rho']);
+ gdat_data.(sub_eff_out).fit.error_bar = cxrs_profiles.proffit.(['d' sub_eff]);
+ else
+ gdat_data.(sub_eff_out).fit.data = [];
+ gdat_data.(sub_eff_out).fit.rho = [];
+ gdat_data.(sub_eff_out).fit.error_bar = [];
+ end
+ % raw data (use all data so keep same size)
+ gdat_data.(sub_eff_out).raw.data = NaN * ones(inb_channels,inb_times);
+ gdat_data.(sub_eff_out).raw.rho = NaN * ones(inb_channels,inb_times);
+ gdat_data.(sub_eff_out).raw.error_bar = NaN * ones(inb_channels,inb_times);
+ gdat_data.(sub_eff_out).raw.error_bar_rho = NaN * ones(inb_channels,inb_times);
+ gdat_data.(sub_eff_out).raw.cxrs_system = NaN * ones(inb_channels,inb_times);
+ gdat_data.time_interval = [];
+ for it=1:inb_times
+ if isfield(cxrs_profiles,sub_eff)
+ nb_raw_points = length(cxrs_profiles.(sub_eff){it}.use.y);
+ gdat_data.(sub_eff_out).raw.data(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.y;
+ gdat_data.(sub_eff_out).raw.rho(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.x;
+ gdat_data.(sub_eff_out).raw.error_bar(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.dy;
+ gdat_data.(sub_eff_out).raw.error_bar_rho(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.dx;
+ gdat_data.(sub_eff_out).raw.cxrs_system(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.sys;
+ gdat_data.time_interval{it} = cxrs_profiles.(sub_eff){it}.t_lim;
+ end
+ end
+ gdat_data.(sub_eff_out).units = sub_nodes_units{i};
+ if i==1
+ gdat_data.data = gdat_data.(sub_eff_out).fit.data;
+ gdat_data.x = gdat_data.(sub_eff_out).fit.rho;
+ gdat_data.error_bar = gdat_data.(sub_eff_out).fit.error_bar;
+ gdat_data.units = gdat_data.(sub_eff_out).units;
+ end
+ end
+ gdat_data.t = cxrs_profiles.proffit.time;
+ gdat_data.dim = {gdat_data.x; gdat_data.t};
+ if isempty(time_interval)
+ gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[1 2 3],[],cxrs_params); % with cxrs_params'];
+ else
+ gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[1 2 3],[' num2str(time_interval) '],cxrs_params); % with cxrs_params'];
+ end
+ gdat_data.dimunits{1} = '';
+ gdat_data.dimunits{2} = 's';
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ 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
+ gdat_data.gdat_params.time = time;
+ if (gdat_params.nverbose>=3); disp(['"time" is expected as an option, choose default time = ' num2str(time)]); end
+ 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;
+ else
+ gdat_data.gdat_params.zshift = 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);
+ if isempty(fnames_readresults)
+ if gdat_params.nverbose>=1;
+ warning(['could not create eqdsk file from read_results_for_chease with data_request= ' data_request_eff]);
+ end
+ return
+ end
+ 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_jet 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;
+ else
+ gdat_data.gdat_params.cocos = cocos_out;
+ 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
+ xx=repmat(reshape(gdat_data.dim{1},length(gdat_data.dim{1}),1),1,size(gdat_data.eqdsk{itime}.psi,2));
+ yy=repmat(reshape(gdat_data.dim{2},1,length(gdat_data.dim{2})),size(gdat_data.eqdsk{itime}.psi,1),1);
+ aa = interpos2Dcartesian(gdat_data.eqdsk{itime}.rmesh,gdat_data.eqdsk{itime}.zmesh ...
+ ,gdat_data.eqdsk{itime}.psi,xx,yy,-1,-1);
+ gdat_data.data(:,:,itime) = aa;
+ end
+ else
+ gdat_data.eqdsk = write_eqdsk(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
+ if shot< 50926
+ n1=tdi('abs(mhdmode("LFS",1,1))');
+ n2=tdi('abs(mhdmode("LFS",2,1))');
+ n3=tdi('abs(mhdmode("LFS",3,1))');
+ gdat_data.data_fullpath='abs(mhdmode("LFS",n,1));% for 1, 2, 3';
+ else
+ if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source)
+ % gdat_data.gdat_params.source;
+ else
+ gdat_data.gdat_params.source = '23';
+ end
+ if length(gdat_data.gdat_params.source)>=2 && strcmp(gdat_data.gdat_params.source(1:2),'23')
+ aaLFSz23_sect3=tdi('\atlas::DT196_MHD_001:channel_067');
+ aaLFSz23_sect11=tdi('\atlas::DT196_MHD_001:channel_075');
+ n1 = aaLFSz23_sect3;
+ n1.data = aaLFSz23_sect3.data - aaLFSz23_sect11.data;
+ n2 = aaLFSz23_sect3;
+ n2.data = aaLFSz23_sect3.data + aaLFSz23_sect11.data;
+ n3=n1;
+ gdat_data.data_fullpath='\atlas::DT196_MHD_001:channel_067 -+ \atlas::DT196_MHD_001:channel_075 for n=1,2, LFS_sect_3/11, z=23cm';
+ if strcmp(gdat_data.gdat_params.source,'23full')
+ % HFS from sec 3 and 11
+ aaHFSz23_sect3=tdi('\atlas::DT196_MHD_002:channel_018');
+ aaHFSz23_sect11=tdi('\atlas::DT196_MHD_002:channel_022');
+ gdat_data.n1_HFS=aaHFSz23_sect3;
+ gdat_data.n1_HFS.data = gdat_data.n1_HFS.data - aaHFSz23_sect11.data;
+ gdat_data.n2_HFS=aaHFSz23_sect3;
+ gdat_data.n2_HFS.data = gdat_data.n2_HFS.data + aaHFSz23_sect11.data;
+ gdat_data.data_fullpath=['\atlas::DT196_MHD_001:channel_067 -+ \atlas::DT196_MHD_001:channel_075 for n=1,2, LFS_sect_3/11, z=23cm; _HFS' ...
+ ' same for sector HFS: MHD_002:channel_018-+MHD_002:channel_022'];
+ end
+ elseif strcmp(gdat_data.gdat_params.source(1),'0')
+ aaLFSz0_sect3=tdi('\atlas::DT196_MHD_001:channel_083');
+ aaLFSz0_sect11=tdi('\atlas::DT196_MHD_001:channel_091');
+ n1 = aaLFSz0_sect3;
+ n1.data = aaLFSz0_sect3.data - aaLFSz0_sect11.data;
+ n2 = aaLFSz0_sect3;
+ n2.data = aaLFSz0_sect3.data + aaLFSz0_sect11.data;
+ n3=n1;
+ gdat_data.data_fullpath='\atlas::DT196_MHD_001:channel_083 -+ \atlas::DT196_MHD_001:channel_091 for n=1,2, LFS_sect_3/11, z=0cm';
+ if strcmp(gdat_data.gdat_params.source,'0full')
+ % sect 11 180deg from sec 3
+ aaHFSz0_sect3=tdi('\atlas::DT196_MHD_002:channel_034');
+ aaHFSz0_sect11=tdi('\atlas::DT196_MHD_002:channel_038');
+ gdat_data.n1_HFS=aaHFSz0_sect11;
+ gdat_data.n1_HFS.data = gdat_data.n1_HFS.data - aaHFSz0_sect11.data;
+ gdat_data.n2_HFS=aaHFSz0_sect11;
+ gdat_data.n2_HFS.data = gdat_data.n2_HFS.data + aaHFSz0_sect11.data;
+ gdat_data.data_fullpath=['\atlas::DT196_MHD_001:channel_083 -+ \atlas::DT196_MHD_001:channel_091 for n=1,2, LFS_sect3/11, z=0cm; _HFS' ...
+ ' same for HFS: MHD_002:channel_034-+MHD_002:channel_038'];
+
+ end
+
+
+ else
+ disp('should not be here in ''mhd'', ask O. Sauter')
+ return
+ end
+ end
+ 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.request_description = 'delta_Bdot from magnetic probes to get n=1, 2 and 3';
+ gdat_data.label = {'n=1','n=2','n=3'}; % can be used in legend(gdat_data.label)
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'ne','te'}
+ % ne or Te from Thomson data on raw z mesh vs (z,t)
+ if ~(isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
+ gdat_data.gdat_params.edge>0)
+ gdat_data.gdat_params.edge = 0;
+ end
+ [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.edge,gdat_params.nverbose);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'ne_rho', 'te_rho', 'nete_rho'}
+ % if nete_rho, do first ne, then Te later (so fir stuff already done)
+ 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
+ if ~(isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
+ gdat_data.gdat_params.edge>0)
+ gdat_data.gdat_params.edge = 0;
+ end
+ [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.edge,gdat_params.nverbose);
+ % construct rho mesh
+ edge_str_ = '';
+ edge_str_dot = '';
+ if gdat_data.gdat_params.edge
+ edge_str_ = '_edge';
+ edge_str_dot = '.edge';
+ end
+ psi_max=gdat([],['\results::thomson' edge_str_dot ':psi_max' substr_liuqe]);
+ psiscatvol=gdat([],['\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
+ if (gdat_params.nverbose>=1);
+ disp(' bad time dimension for zshift')
+ disp(['it should be 1 or ' num2str(length(psiscatvol.dim{1})) ' or ' num2str(length(psitdi.dim{3}))])
+ end
+ 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
+ if gdat_params.nverbose>=1; warning(['psiscatvol empty?, no rho calculated for data_request = ' data_request_eff]); end
+ rho=[];
+ end
+ gdat_data.dim{1}=rho;
+ gdat_data.dimunits=[{'sqrt(psi_norm)'} ; {'time [s]'}];
+ gdat_data.x=rho;
+ %%%%%%%%%%%
+ % 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')
+ % note, now has ne.data_raw for density without fir_to_thomson_ratio
+ gdat_data.ne.data = gdat_data.data;
+ gdat_data.ne.data_raw = gdat_data.data_raw;
+ gdat_data.ne.error_bar = gdat_data.error_bar;
+ gdat_data.ne.error_bar_raw = gdat_data.error_bar_raw;
+ gdat_data.ne.firrat=gdat_data.firrat;
+ gdat_data.ne.units = 'm^{-3}';
+ gdat_data = rmfield(gdat_data,{'firrat','data_raw','error_bar_raw'});
+ %
+ nodenameeff=['\results::thomson' edge_str_dot ':te'];
+ tracetdi=tdi(nodenameeff);
+ nodenameeff=['\results::thomson' edge_str_dot ':te; error_bar'];
+ tracestd=tdi(['\results::thomson' edge_str_dot ':te:error_bar']);
+ gdat_data.te.data=tracetdi.data';
+ gdat_data.te.error_bar=tracestd.data';
+ gdat_data.te.units = tracetdi.units;
+ gdat_data.data_fullpath=['pe=1.6e-19*ne*Te in data, .ne, .te from \results::thomson' ...
+ edge_str_dot ':ne and te and projected on rhopol\_norm'];
+ 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
+ %%%%%%%%%%% add fitted profiles if 'fit'>=1
+ default_fit_type = 'conf';
+ if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit)
+ gdat_data.gdat_params.fit = 1;
+ end
+ % add empty structure for fit so results is always the same with or without call to fit=1 or 0
+ gdat_data.fit.data = [];
+ gdat_data.fit.x = [];
+ gdat_data.fit.t = [];
+ gdat_data.fit.units = [];
+ gdat_data.fit.data_fullpath = [];
+ if strcmp(data_request_eff(1:4),'nete')
+ % note gdat_data.fit.data level is for pe
+ gdat_data.fit.ne.data = [];
+ gdat_data.fit.ne.x = [];
+ gdat_data.fit.ne.t = [];
+ gdat_data.fit.ne.units = [];
+ gdat_data.fit.ne.data_fullpath = [];
+ gdat_data.fit.te.data = [];
+ gdat_data.fit.te.x = [];
+ gdat_data.fit.te.t = [];
+ gdat_data.fit.te.units = [];
+ gdat_data.fit.te.data_fullpath = [];
+ end
+ if gdat_data.gdat_params.fit>0
+ % default is from proffit:avg_time
+ if ~(isfield(gdat_data.gdat_params,'fit_type') && ~isempty(gdat_data.gdat_params.fit_type)) || ~any(strcmp(gdat_data.gdat_params.fit_type,{'local', 'avg', 'conf'}))
+ gdat_data.gdat_params.fit_type = default_fit_type;
+ end
+ switch gdat_data.gdat_params.fit_type
+ case 'avg'
+ def_proffit = '\results::proffit.avg_time:';
+ case 'local'
+ def_proffit = '\results::proffit.local_time:';
+ case 'conf'
+ def_proffit = '\results::conf:';
+ otherwise
+ if (gdat_params.nverbose>=1);
+ disp('should not be in switch gdat_data.gdat_params.fit_type')
+ disp('ask olivier.sauter@epfl.ch')
+ end
+ return
+ end
+ if strcmp(gdat_data.gdat_params.fit_type,'conf')
+ nodenameeff = [def_proffit data_request_eff(1:2)];
+ else
+ 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
+ if (gdat_params.nverbose>=1);
+ disp(['should not be here: data_request_eff, data_request_eff(1:2)= ',data_request_eff, data_request_eff(1:2)]);
+ end
+ end
+ 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)
+ if ~isempty(tracetdi.data) && ~isempty(tracetdi.dim) && ~ischar(tracetdi.data)
+ gdat_data.fit.data = tracetdi.data;
+ else
+ if gdat_params.nverbose>=1
+ disp([nodenameeff ' is empty, thus no fits, check hldsi(shot) and may need to run anaprofs?'])
+ end
+ if strcmp(data_request_eff(1:4),'nete')
+ gdat_data.fit.ne.data_fullpath = [nodenameeff ' is empty'];
+ gdat_data.fit.te.data_fullpath = [nodenameeff ' is empty'];
+ else
+ gdat_data.fit.data_fullpath = [nodenameeff ' is empty'];
+ end
+ return
+ end
+ else
+ if ~isempty(tracetdi.data) && size(tracetdi.data,3)>=trialindx+1
+ gdat_data.fit.data = tracetdi.data(:,:,trialindx+1);
+ else
+ if gdat_params.nverbose>=1
+ disp([nodenameeff ' with trialindx=' num2str(trialindx) ' is empty, thus no fits, check hldsi(shot) and may need to run anaprofs?'])
+ end
+ 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_jet.timedim~=2 | mapping_for_jet.gdat_timedim~=2
+ if (gdat_params.nverbose>=1);
+ disp(['unexpected timedim in fit: data_request_eff= ' data_request_eff ...
+ ', mapping_for_jet.timedim= ' mapping_for_jet.timedim ...
+ ', mapping_for_jet.gdat_timedim= ' mapping_for_jet.gdat_timedim]);
+ end
+ end
+ 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.x = gdat_data.fit.x;
+ gdat_data.fit.ne.t = gdat_data.fit.t;
+ gdat_data.fit.ne.units = gdat_data.fit.units;
+ gdat_data.fit.ne.data_fullpath = gdat_data.fit.data_fullpath;
+ if strcmp(gdat_data.gdat_params.fit_type,'conf')
+ nodenameeff = [def_proffit 'te'];
+ else
+ nodenameeff = [def_proffit 'teft'];
+ end
+ 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
+ gdat_data.fit.te.x = gdat_data.fit.ne.x;
+ gdat_data.fit.te.t = gdat_data.fit.ne.t;
+ if any(strcmp(fieldnames(tracetdi),'units'))
+ gdat_data.fit.te.units=tracetdi.units;
+ end
+ gdat_data.fit.te.data_fullpath = [nodenameeff];
+ % 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_type = default_fit_type;
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'powers'}
+ % note: same time array for all main, ec, ohm, nbi, ...
+ % At this stage fill just ech, later add nbi
+
+ % EC
+ nodenameeff='\results::toray.input:p_gyro';
+ tracetdi=tdi(nodenameeff);
+ if isempty(tracetdi.data) || isempty(tracetdi.dim) % || ischar(tracetdi.data) (to add?)
+ if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
+ gdat_data.ec.data = [];
+ gdat_data.ec.units = [];
+ gdat_data.ec.dim=[];
+ gdat_data.ec.dimunits=[];
+ gdat_data.ec.t=[];
+ gdat_data.ec.x=[];
+ gdat_data.ec.data_fullpath=[];
+ gdat_data.ec.label='';
+ else
+ 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';
+ end
+
+ % Ohmic
+ gdat_data.ohm.data = [];
+ gdat_data.ohm.units = '';
+ gdat_data.ohm.dim=gdat_data.dim;
+ gdat_data.ohm.dimunits=gdat_data.dimunits;
+ gdat_data.ohm.t=[];
+ gdat_data.ohm.x=[];
+ gdat_data.ohm.data_fullpath=[];
+ gdat_data.ohm.label='';
+ % get ohmic power simply from vloop*Ip (minus sign for JET)
+ ip=gdat([],'ip');
+ vloop=gdat([],'vloop');
+ tension = -1e5;
+ if isempty(gdat_data.t)
+ gdat_data.t = vloop.t;
+ gdat_data.dim{1} = gdat_data.t;
+ gdat_data.ohm.data = vloop.data;
+ gdat_data.dimunits{1} = 's';
+ gdat_data.units = 'W';
+ else
+ 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;
+ end
+ gdat_data.ohm.units = gdat_data.units;
+ 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;
+ if ~isempty(gdat_data.ec.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}';
+ else
+ gdat_data.dim{2} = [1];
+ gdat_data.dimunits{2} = 'Pohm=Ptot';
+ gdat_data.data_fullpath=['tot power from ohm'];
+ gdat_data.label = 'P_{tot}=P_{ohm}';
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ 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);
+ if isempty(tracetdi.data) || isempty(tracetdi.dim) % || ischar(tracetdi.data) (to add?)
+ if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
+ if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
+ return
+ end
+ 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);
+ if isempty(tracetdi.data) || isempty(tracetdi.dim) || ~any(~isnan(tracetdi.data)) % || ischar(tracetdi.data) (to add?)
+ if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
+ if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
+ return
+ end
+ 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_jet([],params_eff);
+ if isempty(q_rho.data) || isempty(q_rho.dim) % || ischar(q_rho.data) (to add?)
+ if (gdat_params.nverbose>=1); warning(['problems loading data for q_rho for data_request= ' data_request_eff]); end
+ if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
+ return
+ end
+ params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE
+ psi_axis=gdat_jet([],params_eff);
+ params_eff.data_request='b0'; % psi_edge=0 with LIUQE
+ b0=gdat_jet([],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)
+ if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
+ 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)
+ % try to run psitbxput
+ psitbxput_version = 1.3;
+ psitbxput(psitbxput_version,shot);
+ ishot = mdsopen(shot);
+ tracetdi=tdi(nodenameeff);
+ if isempty(tracetdi.data) || isempty(tracetdi.dim)
+ return
+ end
+ 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
+ if (gdat_params.nverbose>=1)
+ disp(['should not be here in vol cases with data_request = ' data_request_eff]);
+ end
+ return
+ end
+ end
+ gdat_data.t = gdat_data.dim{mapping_for_jet.gdat_timedim};
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'sxr', 'mpx'}
+
+ if strcmp(data_request_eff,'mpx')
+ data_request_eff = 'mpx'; % mpx chosen through parameter 'source' within 'sxr'
+ gdat_data.data_request = data_request_eff;
+ gdat_data.gdat_params.source = 'mpx';
+ end
+ % sxr from dmpx by default or xtomo if 'camera','xtomo' is provided
+ if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
+ gdat_data.gdat_params.source = 'mpx';
+ elseif ~strcmp(lower(gdat_data.gdat_params.source),'xtomo') && ~strcmp(lower(gdat_data.gdat_params.source),'mpx')
+ if gdat_data.gdat_params.nverbose>=1
+ warning(['source = ' gdat_data.gdat_params.source ' not expected with data_request= ' data_request_eff])
+ end
+ return
+ end
+ if ~isfield(gdat_data.gdat_params,'time_interval')
+ gdat_data.gdat_params.time_interval = [];
+ end
+ if ~isfield(gdat_data.gdat_params,'freq')
+ gdat_data.gdat_params.freq = 'slow';
+ end
+ switch lower(gdat_data.gdat_params.source)
+ case {'mpx', 'dmpx'}
+ gdat_data.gdat_params.source = 'mpx';
+ if ~isfield(gdat_data.gdat_params,'camera') || isempty(gdat_data.gdat_params.camera)
+ gdat_data.gdat_params.camera = 'top';
+ else
+ gdat_data.gdat_params.camera = lower(gdat_data.gdat_params.camera);
+ end
+ if ~any(liuqe_version==[1, 2, 3])
+ if gdat_data.gdat_params.nverbose>=3
+ disp(['liuqe_version = ' liuqe_version ' not supported for data_request= ' data_request_eff]);
+ end
+ return
+ end
+ freq_opt = 0;
+ if strcmp(gdat_data.gdat_params.freq,'fast'); freq_opt = 1; end
+ t_int = [0 10]; % starts from 0 otherwise mpxdata gives data from t<0
+ if ~isempty(gdat_data.gdat_params.time_interval); t_int = gdat_data.gdat_params.time_interval; end
+ gdat_data.top.data = [];
+ gdat_data.top.x = [];
+ gdat_data.top.channel = [];
+ gdat_data.bottom.data = [];
+ gdat_data.bottom.x = [];
+ gdat_data.bottom.channel = [];
+ try
+ mpx = mpxdata(shot,'svgr','freq',freq_opt,'liuqe',liuqe_version,'detec',gdat_data.gdat_params.camera, ...
+ 'time',t_int);
+ catch
+ if gdat_data.gdat_params.nverbose>=1
+ warning('problem with mpxdata')
+ end
+ return
+ end
+ gdat_data.units = {'au'}; % not known at this stage
+ gdat_data.dimunits = {'', 's'};
+ gdat_data.data_fullpath = ['using mpxdata(' num2str(shot) ',''svgr'') with params in gdat_data.gdat_params'];
+ if ~strcmp(gdat_data.gdat_params.camera,'both')
+ % invert index of time and channel (rho)
+ gdat_data.data = mpx.(gdat_data.gdat_params.camera(1:3)).signal.data';
+ gdat_data.t = mpx.(gdat_data.gdat_params.camera(1:3)).signal.dim{1};
+ gdat_data.dim{1} = interp1(mpx.top.rho.time,mpx.top.rho.rhopsi,gdat_data.t)';
+ gdat_data.dim{2} = mpx.(gdat_data.gdat_params.camera(1:3)).signal.dim{1};
+ gdat_data.x = gdat_data.dim{1};
+ gdat_data.(gdat_data.gdat_params.camera).data = gdat_data.data;
+ gdat_data.(gdat_data.gdat_params.camera).x = gdat_data.x;
+ gdat_data.(gdat_data.gdat_params.camera).channel = mpx.(gdat_data.gdat_params.camera(1:3)).signal.dim{2};
+ gdat_data.data_fullpath = ['MPX for ' gdat_data.gdat_params.camera ' camera in .data, "rho" in .x between [-1,1]' ...
+ char(10) gdat_data.data_fullpath];
+ else
+ gdat_data.data = mpx.signal.data';
+ gdat_data.t = mpx.signal.dim{1};
+ gdat_data.dim{1} = interp1(mpx.top.rho.time,mpx.top.rho.rhopsi,gdat_data.t);
+ gdat_data.dim{2} = mpx.(gdat_data.gdat_params.camera(1:3)).signal.dim{1};
+ %
+ gdat_data.top.data = mpx.top.signal.data';
+ gdat_data.top.x = gdat_data.dim{1};
+ gdat_data.top.channel = mpx.top.signal.dim{2};
+ gdat_data.bottom.data = mpx.bot.signal.data;
+ gdat_data.bottom.x = interp1(mpx.bot.rho.time,mpx.bot.rho.rhopsi,gdat_data.t);
+ gdat_data.bottom.channel = mpx.bottom.signal.dim{2};
+ gdat_data.(gdat_data.gdat_params.camera).channel = mpx.(gdat_data.gdat_params.camera(1:3)).signal.dim{2};
+ gdat_data.data_fullpath = ['MPX for ' gdat_data.gdat_params.camera ' camera in .data, "rho" in .x between [-1,1]' ...
+ char(10) gdat_data.data_fullpath];
+ gdat_data.x = gdat_data.dim{1};
+ end
+
+ case 'xtomo'
+ % so far allow string and array as 'camera' choices:
+ % camera = [] (default, thus get XTOMOGetData defaults), 'central', [3 6] (camera numbers)
+ camera_xtomo = [];
+ channel_xtomo = [];
+ if ~isfield(gdat_data.gdat_params,'camera') || isempty(gdat_data.gdat_params.camera)
+ gdat_data.gdat_params.camera = [];
+ elseif isnumeric(gdat_data.gdat_params.camera)
+ if length(gdat_data.gdat_params.camera) > 10
+ if gdat_data.gdat_params.nverbose>=1; warning('max number of camera is 10 for XTOMO'); end
+ gdat_data.gdat_params.camera = gdat_data.gdat_params.camera(1:10);
+ end
+ camera_xtomo = zeros(1,10);
+ camera_xtomo(gdat_data.gdat_params.camera) = 1;
+ elseif ischar(gdat_data.gdat_params.camera)
+ gdat_data.gdat_params.camera = lower(gdat_data.gdat_params.camera);
+ if strcmp(gdat_data.gdat_params.camera,'central')
+ camera_xtomo = zeros(1,10);
+ icam = 3
+ camera_xtomo(icam) = 1;
+ channel_xtomo = (icam-1)*20 + 9;
+ end
+ else
+ if gdat_data.gdat_params.nverbose>=3; disp(['camera = ' gdat_data.gdat_params.camera ' not implemented']); end
+ gdat_data.gdat_params.camera = [];
+ end
+
+ try
+ if isempty(gdat_data.gdat_params.time_interval);
+ [sig,t,Sat_channel,cat_default] = XTOMOGetData(shot,0,4,camera_xtomo,[]);
+ else
+ [sig,t,Sat_channel,cat_default] = XTOMOGetData(shot,gdat_data.gdat_params.time_interval(1), ...
+ gdat_data.gdat_params.time_interval(2),camera_xtomo,[]);
+ end
+ catch
+ if gdat_data.gdat_params.nverbose>=1
+ warning('problem with XTOMOGetData, no data')
+ end
+ return
+ end
+ gdat_data.t = t;
+ gdat_data.units = 'au';
+ gdat_data.xtomo.extra_params = cat_default;
+ gdat_data.xtomo.saturated_channels = Sat_channel;
+ gdat_data.data_fullpath = ['using XTOMOGetData(' num2str(shot) ') with some additional params from gdat_data.gdat_params'];
+ if isempty(channel_xtomo)
+ % provide all chords
+ gdat_data.data = sig;
+ gdat_data.x = [1:size(sig,1)];
+ gdat_data.dimunits = {'20 chords per camera'; 's'};
+ else
+ keyboard
+ % extract only given channels
+ gdat_data.data = sig(channel_xtomo,:);
+ gdat_data.x = channel_xtomo;
+ gdat_data.dimunits = {'chords where floor(dim{1}/10) gives camera nb and mod(dim{1},10) chord nb'; 's'};
+ end
+ gdat_data.dim = {gdat_data.x; gdat_data.t};
+
+ otherwise
+ if gdat_data.gdat_params.nverbose>=1
+ warning(['camera = ' gdat_data.gdat_params.camera ' not expected with data_request= ' data_request_eff])
+ end
+ return
+
+ end
+
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'profnerho','profterho'}
+ % for backward compatibility but corresponds to ne_rho with param.fit_type='auto' (JET 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 gdat_params.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
+ if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_jet']); end
+ error_status=901;
+ return
+ end
+
+else
+ if (gdat_params.nverbose>=1); warning(['JET method=' mapping_for_jet.method ' not known yet, contact Olivier.Sauter@epfl.ch']); end
+ error_status=602;
+ return
+end
+
+%if ishot==shot; mdsclose; end
+
+gdat_data.gdat_params.help = jet_help_parameters(fieldnames(gdat_data.gdat_params));
+gdat_data.mapping_for.jet = mapping_for_jet;
+gdat_params = gdat_data.gdat_params;
+error_status=0;
+
+return
+
diff --git a/crpptbx/JET/jet_help_parameters.m b/crpptbx/JET/jet_help_parameters.m
new file mode 100644
index 0000000000000000000000000000000000000000..2167c5fc34ca10a4e283eb899803589a07df34ef
--- /dev/null
+++ b/crpptbx/JET/jet_help_parameters.m
@@ -0,0 +1,70 @@
+function help_struct = jet_help_parameters(parameter_list)
+%
+% retrieve from present table the relevant help lines for the parameter_list{:}
+% should come from sqlite database at some point...
+%
+% return the whole help structure if parameter_list empty or not provided
+%
+% do:
+% help_struct = jet_help_parameters(fieldnames(gdat_data.gdat_params));
+%
+% to get relevant help description
+%
+
+% Defaults
+help_struct_all = struct(...
+ 'data_request', ['automatically filled in by gdat, name of request used in gdat call.' char(10) ...
+ 'contains current list of keywords if gdat called with no arguments: aa=gdat;' char(10) ...
+ 'Note shot value should not be in params so params can be used to load same data from another shot'] ...
+ ,'machine', 'machine name like ''TCV'', ''JET'', case insensitive' ...
+ ,'doplot', '0 (default), if 1 calls gdat_plot for a new figure, -1 plot over current figure with hold all, see gdat_plot for details' ...
+ ,'liuqe','liuqe version 1 (default), 2, 3 for LIUQE1, 2, 3 resp. or -1 for model values' ...
+ ,'nverbose','1 (default) displays warnings, 0: only errors, >=3: displays all extra information' ...
+ );
+
+% JET related
+% $$$ help_struct_all.cxrs_plot = '0 (default) no plots, 1 get plots from CXRS_get_profiles see ''help CXRS_get_profiles'' for k_plot values';
+% $$$ help_struct_all.time_interval = ['if provided sets a specific time interval [tstart tend].' ...
+% $$$ char(10) 'cxrs: (time_interval can have several nbs) take data and average over time interval(s) only, plots from CXRS_get_profiles are then provided' ...
+% $$$ ' as well'];
+% $$$ help_struct_all.fit_tension = ['smoothing value used in interpos fitting routine, -30 means ''30 times default value'', thus -1 often a' ...
+% $$$ ' good value' char(10) ...
+% $$$ 'cxrs: if numeric, default for all cases, if structure, default for non given fields'];
+help_struct_all.time = 'time(s) value(s) if relevant, for example eqdsk is provided by default only for time=1.0s';
+% $$$ help_struct_all.zshift = 'vertical shift of equilibrium, either for eqdsk (1 to shift to zaxis=0) or for mapping measurements on to rho surfaces [m]';
+help_struct_all.cocos = ['cocos value desired in output, uses eqdsk_cocos_transform. Note should use latter if a specific Ip and/or B0 sign' ...
+ 'is wanted. See O. Sauter et al Comput. Phys. Commun. 184 (2013) 293'];
+% $$$ help_struct_all.edge = '0 (default), 1 to get edge Thomson values';
+% $$$ help_struct_all.fit = '0, no fit profiles, 1 (default) if fit profiles desired as well, relevant for _rho profiles. See also fit_type';
+% $$$ help_struct_all.fit_type = 'provenance of fitted profiles ''conf'' (default) from conf nodes, ''avg'' or ''local'' for resp. proffit: nodes';
+help_struct_all.source = 'sxr: ''G'' (default, with ssx), camera name ''J'', ''G'', ...[[F-M], case insensitive';
+help_struct_all.camera = ['[] (default, all), [i1 i2 ...] chord nbs ([1 3 5] if only chords 1, 3 and 5 are desired), ''central'' uses J_049'];
+help_struct_all.freq = '''slow'', default (means ssx, 500kHz), lower sampling; ''fast'' full samples 2MHz; integer value nn for downsampling every nn''th points';
+%help_struct_all. = '';
+
+%help_struct_all. = '';
+
+
+
+if ~exist('parameter_list') || isempty(parameter_list)
+ help_struct = help_struct_all;
+ return
+end
+
+if iscell(parameter_list)
+ parameter_list_eff = parameter_list;
+elseif ischar(parameter_list)
+ % assume only one parameter provided as char string
+ parameter_list_eff{1} = parameter_list;
+else
+ disp(['unknown type of parameter_list in jet_help_parameters.m'])
+ parameter_list
+ return
+end
+
+fieldnames_all = fieldnames(help_struct_all);
+for i=1:length(parameter_list_eff)
+ if any(strcmp(fieldnames_all,parameter_list_eff{i}))
+ help_struct.(parameter_list_eff{i}) = help_struct_all.(parameter_list_eff{i});
+ end
+end
diff --git a/crpptbx/JET/jet_requests_mapping.m b/crpptbx/JET/jet_requests_mapping.m
new file mode 100644
index 0000000000000000000000000000000000000000..97a3fc86ddf7c171830df2a760916c7e636aa4ae
--- /dev/null
+++ b/crpptbx/JET/jet_requests_mapping.m
@@ -0,0 +1,323 @@
+function mapping = jet_requests_mapping(data_request)
+%
+% Information pre-defined for gdat_jet, you can add cases here to match official cases in gdat_jet, 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 JET, following choices are set so far:
+% method = 'signal' then expression contains the shotfile, diagnostic and if needed the experiment
+% expression is a cell array
+% method = 'expression', then expression is the expression to return gdat_tmp...
+% method = 'switchcase', then there will be a specific case within gdat_jet (usual case when not directly a signal)
+%
+% label is used for plotting
+if iscell(data_request) % || (~ischar(data_request) && length(data_request)>1)
+ mapping.label = data_request;
+ mapping.method = 'signal'; % assume a full tracename is given, so just try with tdi (could check there are some ":", etc...)
+ mapping.expression = data_request;
+ mapping.gdat_timedim = mapping.timedim;
+ return
+end
+
+switch lower(data_request)
+ case 'a_minor'
+ mapping.timedim = 1;
+ mapping.label = 'a\_minor';
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''r_inboard'';' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);gdat_tmp.r_inboard=gdat_tmp.data;' ...
+ 'params_eff.data_request=''r_outboard'';' ...
+ 'gdat_tmp2=gdat_jet(shot,params_eff);gdat_tmp.r_outboard=gdat_tmp2.data;' ...
+ 'gdat_tmp.data = 0.5.*(gdat_tmp2.data-gdat_tmp.data);gdat_tmp.label=''' mapping.label ''';' ...
+ 'gdat_tmp.gdat_request=''' data_request ''';'];
+ case 'b0'
+ mapping.timedim = 1;
+ mapping.label = 'B_0';
+ mapping.method = 'signal';
+ mapping.expression = [{'FPC'},{'BTF'}];
+ case 'beta'
+ mapping.timedim = 1;
+ mapping.label = '\beta';
+ mapping.method = 'signal';
+ mapping.expression = [{'TOT'},{'beta'}];
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''betan'';' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);' ...
+ 'params_eff.data_request=''ip'';gdat_tmp2=gdat_jet(shot,params_eff);' ...
+ 'params_eff.data_request=''b0'';gdat_tmp3=gdat_jet(shot,params_eff);' ...
+ 'params_eff.data_request=''a_minor'';gdat_tmp4=gdat_jet(shot,params_eff);' ...
+ 'tmp_data_ip=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ...
+ 'tmp_data_b0=interp1(gdat_tmp3.t,gdat_tmp3.data,gdat_tmp.t,[],NaN);' ...
+ 'tmp_data_a=interp1(gdat_tmp4.t,gdat_tmp4.data,gdat_tmp.t,[],NaN);' ...
+ 'gdat_tmp.data = 0.01.*abs(gdat_tmp.data.*tmp_data_ip./1e6./tmp_data_a./tmp_data_b0);' ...
+ 'gdat_tmp.label=''' mapping.label ''';gdat_tmp.gdat_request=''' data_request ''';'];
+ case 'betan'
+ mapping.timedim = 1;
+ mapping.label = '\beta_N';
+ mapping.method = 'signal';
+ mapping.expression = [{'TOT'},{'beta_N'}];
+ case 'betap'
+ mapping.timedim = 1;
+ mapping.label = '\beta_p';
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'betpol'}];
+ case 'cxrs'
+ mapping.timedim = 2;
+ mapping.label = 'cxrs';
+ mapping.method = 'switchcase';
+ mapping.expression = '';
+ case 'delta'
+ mapping.timedim = 1;
+ mapping.label = 'delta';
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''delta_bottom''; ' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);params_eff.data_request=''delta_top'';' ...
+ 'gdat_tmp2=gdat_jet(shot,params_eff);gdat_tmp.data = 0.5.*(gdat_tmp.data+gdat_tmp2.data);'];
+ case 'delta_top'
+ mapping.label = 'delta\_top';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'delRoben'}];
+ case 'delta_bottom'
+ mapping.label = 'delta\_bottom';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'delRuntn'}];
+ case {'ece', 'eced', 'ece_rho', 'eced_rho'}
+ mapping.timedim = 2;
+ mapping.method = 'switchcase';
+ mapping.expression = '';
+ case 'eqdsk'
+ mapping.timedim = 2;
+ mapping.method = 'switchcase'; % could use function make_eqdsk directly?
+ mapping.expression = '';
+ case 'equil'
+ mapping.gdat_timedim = 2;
+ mapping.method = 'switchcase'; % could use function make_eqdsk directly?
+ mapping.expression = '';
+ case 'halpha'
+ mapping.timedim = 1;
+ mapping.label = 'Halpha';
+ mapping.method = 'signal';
+ mapping.expression = [{'POT'},{'ELMa-Han'}];
+ case 'ioh'
+ mapping.timedim = 1;
+ mapping.label = 'I ohmic transformer';
+ mapping.method = 'signal';
+ mapping.expression = [{'MBI'},{'IOH'}];
+ case 'ip'
+ mapping.timedim = 1;
+ mapping.label = 'Plasma current';
+ mapping.method = 'signal';
+ mapping.expression = [{'ppf'},{'efit'},{'xip'}];
+ mapping.units = 'A';
+ case 'kappa'
+ mapping.timedim = 1;
+ mapping.label = '\kappa';
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'k'}];
+ case 'mhd'
+ mapping.timedim = 1;
+ mapping.label = 'Odd and Even n';
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request={''MOD'',''OddN''}; ' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);gdat_tmp.data=reshape(gdat_tmp.data,length(gdat_tmp.data),1 );' ...
+ 'gdat_tmp.dim{1}=gdat_tmp.t;gdat_tmp.dim{2}=[1 2];gdat_tmp.x=gdat_tmp.dim{2};' ...
+ 'gdat_tmp.n_odd.data = gdat_tmp.data;gdat_tmp.n_odd.data_request=params_eff.data_request;' ...
+ 'params_eff.data_request={''MOD'',''EvenN''};' ...
+ 'gdat_tmp2=gdat_jet(shot,params_eff);gdat_tmp.data(:,2)=reshape(gdat_tmp2.data,length(gdat_tmp2.data),1);' ...
+ 'gdat_tmp.n_even.data = gdat_tmp2.data;gdat_tmp.n_even.data_request=params_eff.data_request;gdat_tmp.label={''n\_odd'',''n\_even''};' ...
+ 'gdat_tmp.full_path=''MOD/Odd in data and .n_odd; .n_even'';' ...
+ 'gdat_tmp.gdat_request=''mhd'';gdat_tmp.gdat_params.data_request=gdat_tmp.gdat_request;'];
+ case 'ne'
+ mapping.timedim = 2;
+ mapping.method = 'switchcase';
+ case 'neint'
+ mapping.timedim = 1;
+ mapping.label = 'line integrated el. density';
+ mapping.method = 'signal';
+ mapping.expression = [{'DCN'},{'H-1'},{'JETD'}];
+ case 'nel'
+ mapping.timedim = 1;
+ mapping.label = 'line-averaged el. density';
+ mapping.expression = [{'FPG'},{'lenH-1'},{'JETD'}];
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''neint'';' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);params_eff.data_request=[{''FPG''},{''lenH-1''},{''JETD''}];' ...
+ 'gdat_tmp2=gdat_jet(shot,params_eff);ij=find(gdat_tmp2.data==0);gdat_tmp2.data(ij)=NaN;' ...
+ 'tmp_data=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ...
+ 'gdat_tmp.data = gdat_tmp.data./(tmp_data+1e-5);'];
+ case 'ne_rho'
+ mapping.timedim = 2;
+ mapping.label = 'ne';
+ mapping.method = 'switchcase';
+ case 'nete_rho'
+ mapping.timedim = 2;
+ mapping.label = 'ne and Te';
+ mapping.method = 'switchcase';
+ case 'ni'
+ 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 = 'switchcase'; % there is psi_axis-psi_edge in FPG but otherwise complicated to get from equil, thus needs swticth case
+ mapping.label ='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.label = 'q_0';
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'q0'},{'JETD'}];
+ case 'q95'
+ mapping.timedim = 1;
+ mapping.label = 'q_{95}';
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'q95'},{'JETD'}];
+ case 'q_edge'
+ mapping.timedim = 1;
+ mapping.label = 'q_{edge}}';
+ mapping.method = 'expression';
+ mapping.method = 'switchcase';
+ mapping.expression = [{'FPG'},{'q95'},{'JETD'}];
+ mapping.expression = [];
+ case 'q_rho'
+ mapping.timedim = 2;
+ mapping.gdat_timedim = 2;
+ mapping.label = 'q';
+ mapping.method = 'switchcase';
+ case 'rgeom'
+ mapping.label = 'Rgeom';
+ mapping.timedim = 1;
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''r_inboard'';' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);gdat_tmp.r_inboard=gdat_tmp.data;' ...
+ 'params_eff.data_request=''r_outboard'';' ...
+ 'gdat_tmp2=gdat_jet(shot,params_eff);gdat_tmp.r_outboard=gdat_tmp2.data;' ...
+ 'gdat_tmp.data = 0.5.*(gdat_tmp2.data+gdat_tmp.data);gdat_tmp.label=''' mapping.label ''';' ...
+ 'gdat_tmp.gdat_request=''' data_request ''';'];
+ case 'r_inboard'
+ mapping.label = 'R\_inboard';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'Rin'},{'JETD'}];
+ case 'r_outboard'
+ mapping.label = 'R\_outboard';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'Raus'},{'JETD'}];
+ case 'rhotor'
+ mapping.timedim = 2;
+ mapping.method = 'switchcase';
+ mapping.label = 'rhotor\_norm';
+ case 'rhotor_edge'
+ mapping.timedim = 1;
+ mapping.method = 'switchcase';
+ mapping.label = 'rhotor\_edge';
+ case 'rhovol'
+ mapping.timedim = 2;
+ mapping.label = 'rhovol\_norm';
+ mapping.method = 'switchcase';
+ case 'rmag'
+ mapping.label = 'R\_magaxis';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'Rmag'},{'JETD'}];
+ case 'sxr'
+ mapping.timedim = 1;
+ mapping.gdat_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 'ti'
+ mapping.label = 'Ti';
+ mapping.method = 'switchcase';
+ case 'vloop'
+ mapping.label = 'Vloop';
+ mapping.timedim = 1;
+ % mapping.method = 'signal';
+ % mapping.expression = [{'MAG'},{'ULid12'},{'JETD'}];
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=[{''MAG''},{''ULid12''},{''JETD''}];' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);ij=find(~isnan(gdat_tmp.data));' ...
+ 'tmp_data=interpos(gdat_tmp.t,gdat_tmp.data,-3e4);' ...
+ 'gdat_tmp.data_smooth = tmp_data;gdat_tmp.gdat_request=''vloop'';gdat_tmp.gdat_params.data_request=gdat_tmp.gdat_request;'];
+ case 'volume'
+ mapping.label = 'Volume';
+ mapping.timedim = 1;
+ mapping.method = 'switchcase';
+ mapping.expression = [{'FPG'},{'Vol'},{'JETD'}];
+ case 'volume_rho'
+ mapping.timedim = 2;
+ mapping.label = 'volume\_norm';
+ mapping.method = 'switchcase';
+ case 'zeff'
+ mapping.label = 'zeff from cxrs';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'ZES'},{'Zeff'},{'JETD'}];
+ case 'zgeom'
+ mapping.label = 'Zgeom';
+ mapping.timedim = 1;
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=[{''FPG''},{''Zoben''},{''JETD''}];' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);gdat_tmp.z_top=gdat_tmp.data;' ...
+ 'params_eff.data_request=[{''FPG''},{''Zunt''},{''JETD''}];' ...
+ 'gdat_tmp2=gdat_jet(shot,params_eff);gdat_tmp.z_bottom=gdat_tmp2.data;' ...
+ 'gdat_tmp.data = 0.5.*(gdat_tmp2.data+gdat_tmp.data);gdat_tmp.label=''' mapping.label ''';' ...
+ 'gdat_tmp.gdat_request=''' data_request ''';'];
+
+ case 'zmag'
+ mapping.label = 'Z\_magaxis';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'FPG'},{'Zmag'},{'JETD'}];
+ %
+ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % extra JET cases (not necessarily in official data_request name list)
+ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ case 'transp'
+ mapping.label = 'transp output';
+ mapping.method = 'switchcase';
+
+
+ otherwise
+ mapping.label = data_request;
+ mapping.method = 'signal'; % assume a full tracename is given, so just try with tdi (could check there are some ":", etc...)
+ mapping.expression = data_request;
+
+end
+
+if isempty(mapping.gdat_timedim)
+ mapping.gdat_timedim = mapping.timedim;
+end
diff --git a/crpptbx/JET/loadJETdata.m b/crpptbx/JET/loadJETdata.m
index c2e4ab03001810dbc08899344e3df3fde564e7e3..7984f23cfe91e298d4d77f4c04ffad25ce341dc5 100644
--- a/crpptbx/JET/loadJETdata.m
+++ b/crpptbx/JET/loadJETdata.m
@@ -321,6 +321,7 @@ switch JETkeywrdcase{index}
end
ij=find(tracename~='''');
tracename=tracename(ij);
+keyboard
[a,x,t,d,e]=rda_eff(shot,ppftype,tracename);
switch tracename
case {'efit/btpd','efit/btpd?uid=jetppf+seq=0'}
diff --git a/crpptbx/JET/rda_jet.m b/crpptbx/JET/rda_jet.m
new file mode 100644
index 0000000000000000000000000000000000000000..8b00c2a2fdbf0cbcab987e5fd8de88bb61859cbd
--- /dev/null
+++ b/crpptbx/JET/rda_jet.m
@@ -0,0 +1,189 @@
+function [data,x,time,hsig,error]=rda_jet(shot,pftype,tracename,varargin);
+%
+% gets data using RDA or mdsplus
+% 1D arrays: assumes dimension is time
+% 2D arrays: assumes data vs (x,time)
+% 3D arrays: assumes data vs (x,time,hsig) (for mdsplus)
+%
+% varargin{1}: time interval or timevalue, will get data closest to that time or within that time interval
+% (DOES NOT WORK YET)
+%
+% examples:
+% [data,x,time,hsig,error]=rda_jet(51994,'ppf','efit/xip');
+% [data,x,time,hsig,error]=rda_jet(52206,'ppf','equi/rmji?uid=jetthg+seq=122');
+%
+% set global variable: usemdsplus to decide if RDA or mdsplus is used:
+% >> global usemdsplus
+% >> usemdsplus=1 % means use mds to get data (default if not defined)
+% >> usemdsplus=0 % means use jetreaddata routine (RDA)
+% if ~exist('usemdsplus'); usemdsplus=1; end
+%
+
+global usemdsplus
+if isempty(usemdsplus); usemdsplus=1; end
+time_int=[];
+if nargin>=4 & ~isempty(varargin{1})
+ time_int=varargin{1};
+end
+
+if usemdsplus
+ % use mdsplus
+
+ if ~unix('test -d /home/duval/mdsplus')
+ addpath('/home/duval/mdsplus')
+ end
+ mdsconnect('mdsplus.jet.efda.org');
+
+ % defines trace to fetch
+ % after '?' specific details
+ separator='+';
+ mainseparator='?';
+ imaintrace=findstr(mainseparator,tracename);
+ if isempty(imaintrace)
+ maintrace=tracename;
+ uid=[];
+ seq=[];
+ diag=[];
+ type=[];
+ else
+ maintrace=tracename(1:imaintrace-1);
+ rest=tracename(imaintrace+1:end);
+ % gets uid if any
+ iuid=findstr('uid=',rest);
+ if isempty(iuid)
+ uid=[];
+ else
+ ii=findstr(separator,rest(iuid:end));
+ if isempty(ii)
+ uid=rest(iuid+4:end);
+ else
+ uid=rest(iuid+4:iuid+ii(1)-2);
+ end
+ end
+ % gets seq if any
+ iseq=findstr('seq=',rest);
+ if isempty(iseq)
+ seq=[];
+ else
+ ii=findstr(separator,rest(iseq:end));
+ if isempty(ii)
+ seq=rest(iseq+4:end);
+ else
+ seq=rest(iseq+4:iseq+ii(1)-2);
+ end
+ end
+ % gets type if any
+ itype=findstr('type=',rest);
+ if isempty(itype)
+ type=[];
+ else
+ ii=findstr(separator,rest(itype:end));
+ if isempty(ii)
+ type=rest(itype+5:end);
+ else
+ type=rest(itype+5:itype+ii(1)-2);
+ end
+ end
+ % gets diag if any
+ idiag=findstr('diag=',rest);
+ if isempty(idiag)
+ diag=[];
+ else
+ ii=findstr(separator,rest(idiag:end));
+ if isempty(ii)
+ diag=rest(idiag+5:end);
+ else
+ diag=rest(idiag+5:idiag+ii(1)-2);
+ end
+ end
+
+ end
+
+ % fetch value
+ if ~isempty(uid)
+ eval(['u=mdsvalue(''_sig=ppfuid("' uid '")'');'])
+ end
+ if strcmpi(type,'lpf')
+ pftype=[type '/' diag];
+ end
+ traceeff=[pftype '/' maintrace];
+ if ~isempty(seq)
+ traceeff=[traceeff '/' num2str(seq)];
+ end
+ user=getenv('USER');
+ eval(['[data,error]=mdsvalue(''_rdaeff' user '=jet("' traceeff '",' num2str(shot) ')'');'])
+ hsig=[];
+ ss=size(data);
+ nbofdim=length(ss);
+ if ss(end)==1; nbofdim=nbofdim-1; end
+ nbofdim=max(nbofdim,1);
+ switch nbofdim
+ case 1
+ eval(['time=mdsvalue(''dim_of(_rdaeff' user ',0)'');']);
+ x=[];
+ if isempty(time) & length(data)>1e6 & strcmpi(type,'lpf') & strcmpi(diag,'kc1f')
+ mdsdisconnect;
+ mdsconnect('mdsplus.jet.efda.org');
+ eval(['aaa=mdsvalue(''_tc91=jet("jpf/da/c1-tc91",' num2str(shot) ');1'');'])
+ taaa=mdsvalue('_ttc91=dim_of(_tc91,0);_ttc91[0]');
+ time=linspace(taaa+1e-6,taaa+4,length(data))';
+ end
+ if isempty(time) & length(data)>1e6 & strcmpi(type,'lpf') & strcmpi(diag,'cats1')
+ ichannel=findstr(':00',maintrace);
+ iblock=str2num(maintrace(ichannel+3));
+ mdsdisconnect;
+ mdsconnect('mdsplus.jet.efda.org');
+ taaa=39.9989+(iblock-1)*8;
+ time=linspace(taaa,taaa+8-4e-6,length(data))';
+ end
+ case 2
+ eval(['x=mdsvalue(''dim_of(_rdaeff' user ',0)'');']);
+ eval(['time=mdsvalue(''dim_of(_rdaeff' user ',1)'');']);
+
+ case 3
+ eval(['x=mdsvalue(''dim_of(_rdaeff' user ',0)'');']);
+ eval(['time=mdsvalue(''dim_of(_rdaeff' user ',1)'');']);
+ disp('3rd dimension in hsig!!!!!!!!!!!!!!!!!!!!!!!!!')
+ eval(['hsig=mdsvalue(''dim_of(_rdaeff' user ',2)'');']);
+
+ otherwise
+ disp([' more than 3 dimensions for ' num2str(shot) ' ; ' pftype '/' tracename])
+ error('in rda_jet')
+
+ end
+
+ mdsdisconnect;
+ if ~unix('test -d /home/duval/mdsplus')
+ rmpath('/home/duval/mdsplus')
+ end
+
+else
+ % use RDA
+ [a,time,x,hsig,error]=jetreaddata(['http://data.jet.uk/' pftype '/' num2str(shot) '/' tracename]);
+ % transpose data as output in C format, reversed from Fortran and matlab standard
+ ss=size(a);
+ nbofdim=length(ss);
+ if ss(end)==1; nbofdim=nbofdim-1; end
+ nbofdim=max(nbofdim,1);
+ if nbofdim==1
+ data=a;
+ else
+ data=a';
+ end
+end
+
+% to prevent problems when trace empty and time become string
+if ischar(time)
+ time=[];
+end
+if ischar(x)
+ x=[];
+end
+if isempty(x) & ~isempty(data) & data==0
+ data=[];
+end
+
+if isempty(data)
+ x=[];
+ time=[];
+end