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 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(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 main default parameters structure
gdat_params.data_request = '';
default_machine = 'tcv';
gdat_params.machine=default_machine;
gdat_params.doplot = 0;
gdat_params.liuqe = 1;
gdat_params.nverbose = 1;
% construct list of keywords from global set of keywords and specific TCV set
% get data_request names from centralized table
data_request_names = get_data_request_names;
% add TCV specific to all:
if ~isempty(data_request_names.tcv)
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);
% 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
gdat_data.gdat_params.help = tcv_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 = mdsipmex(2,'$SHOT');
gdat_data.shot = shot;
do_mdsopen_mdsclose = 0;
if ischar(shot) || isempty(shot)
if gdat_params.nverbose>=1
if isstruct(data_request) && isfield(data_request,'data_request')
warning(['shot cannot be opened with ' data_request.data_request]);
elseif ischar(data_request)
warning(['shot cannot be opened with ' data_request]);
else
warning(['shot cannot be opened']);
end
end
return
end
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
liuqe_version = 1;
if isfield(gdat_data.gdat_params,'liuqe') && ~isempty(gdat_data.gdat_params.liuqe)
liuqe_version = gdat_data.gdat_params.liuqe;
else
gdat_data.gdat_params.liuqe = liuqe_version;
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 || liuqe_version==-1
% 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 isempty(ishot) || ishot~=shot
if gdat_params.nverbose>=1; warning(['cannot open shot= ' num2str(shot)]); end
return
end
end
% fill again at end to have full case, but here to have present status in case of early return
gdat_data.gdat_params.help = tcv_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_data.mapping_for.tcv = mapping_for_tcv;
gdat_params = gdat_data.gdat_params;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 (gdat_params.nverbose>=1); warning(['problems loading data for ' eval_expr ' for data_request= ' data_request_eff]); end
if (gdat_params.nverbose>=3); disp('check .gdat_request list'); end
return
end
gdat_data.data = aatmp.data;
gdat_data.dim = aatmp.dim;
nbdims = length(gdat_data.dim);
if mapping_for_tcv.timedim==-1;
% try to find time dim from units
idim_non1 = []; len_dim = [];
for i=1:length(aatmp.dimunits)
if strcmp(aatmp.dimunits{i},'s'); mapping_for_tcv.timedim = i; end
if length(aatmp.dim{i})>1
idim_non1(end+1) = i;
len_dim(end+1) = length(aatmp.dim{i});
end
end
if length(idim_non1)==1
% only one dim non 1, assume it is time
mapping_for_tcv.timedim = idim_non1(1);
else
[aamax,iaamax]=max(len_dim);
if aamax./min(len_dim)>100
mapping_for_tcv.timedim = idim_non1(iaamax);
end
end
if mapping_for_tcv.timedim==-1
% assume last one except if of length 1
mapping_for_tcv.timedim = nbdims;
if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_tcv.timedim = nbdims-1; end
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 an expression to evaluate and which creates a local structure into variable gdat_tmp
% we copy the structure, to make sure default nodes are defined and to avoid if return is an closed object like tdi
% eval_expr = [mapping_for_tcv.expression ';'];
eval([mapping_for_tcv.expression ';']);
if isempty(gdat_tmp) || (~isstruct(gdat_tmp) & ~isobject(gdat_tmp))
if (gdat_params.nverbose>=1); warning(['expression does not create a gdat_tmp structure: ' mapping_for_tcv.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_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 "expression"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
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_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?)
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_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)
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_tcv';
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_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;
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_tcv.timedim~=2 | mapping_for_tcv.gdat_timedim~=2
if (gdat_params.nverbose>=1);
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
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 TCV)
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_tcv([],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_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)
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_tcv.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' (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 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_tcv']); end
error_status=901;
return
end
else
if (gdat_params.nverbose>=1); warning(['TCV method=' mapping_for_tcv.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 = tcv_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_data.mapping_for.tcv = mapping_for_tcv;
gdat_params = gdat_data.gdat_params;
error_status=0;
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [firthomratio] = get_fir_thom_rat_data(shot,maintracename,timebase,nverbose);
%
% 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)
if (nverbose>=1)
disp('need a timebase in get_fir_thom_rat_data')
end
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)
if (nverbose>=1); disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat'); end
tracefirrat=tdi('\results::thomson:fir_thom_rat');
end
else
tracefirrat=tdi('\results::thomson:fir_thom_rat');
if isempty(tracefirrat.data) || ischar(tracefirrat.data)
if (nverbose>=1); disp('problem with \results::thomson:fir_thom_rat: empty'); end
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)
if (nverbose>=1); disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat'); end
tracefirrat=tdi('\results::thomson:fir_thom_rat');
end
else
tracefirrat=tdi('\results::thomson_edge:fir_thom_rat');
if isempty(tracefirrat.data) || ischar(tracefirrat.data)
if (nverbose>=1); disp('problem with \results::thomson_edge:fir_thom_rat: empty'); end
end
end
else
if (nverbose>=1); disp('bad input in get_fir_thom_rat_data'); end
return
end
if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) && any(~isnan(tracefirrat.data)) ...
&& ~isempty(tracefirrat.dim) && ~isempty(tracefirrat.dim{1})
firthomratio = interp1(tracefirrat.dim{1},tracefirrat.data,timebase);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,doedge,nverbose);
%
try
time=mdsdata('\results::thomson:times');
catch
gdat_data.error_bar = [];
if strcmp(data_request_eff(1:2),'ne')
tracefirrat_data = [];
gdat_data.firrat=tracefirrat_data;
gdat_data.data_raw = gdat_data.data;
gdat_data.error_bar_raw = gdat_data.error_bar;
end
if (nverbose>=1) && shot<100000
warning('Problems with \results::thomson:times')
warning(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
end
return
end
if isempty(time) || ischar(time)
thomsontimes=time;
if (nverbose>=1) && shot<100000
warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check')
disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
end
return
end
edge_str_ = '';
edge_str_dot = '';
if doedge
edge_str_ = '_edge';
edge_str_dot = '.edge';
end
nodenameeff = ['\results::thomson' edge_str_dot ':' data_request_eff(1:2)];
tracetdi=tdi(nodenameeff);
if isempty(tracetdi.data) || ischar(tracetdi.data) || isempty(tracetdi.dim)
gdat_data.error_bar = [];
gdat_data.firrat = [];
gdat_data.data_raw = [];
gdat_data.error_bar_raw = [];
return
end
gdat_data.data=tracetdi.data'; % Thomson data as (t,z)
tracestd=tdi(['\results::thomson' edge_str_dot ':' data_request_eff(1:2) ':error_bar']);
gdat_data.error_bar=tracestd.data';
gdat_data.data_fullpath=[nodenameeff '; error_bar'];
% add fir if ne requested
if strcmp(data_request_eff(1:2),'ne')
tracefirrat_data = get_fir_thom_rat_data(shot,['thomson' edge_str_],time,nverbose);
gdat_data.firrat=tracefirrat_data;
gdat_data.data_raw = gdat_data.data;
gdat_data.data = gdat_data.data_raw * diag(tracefirrat_data);
gdat_data.error_bar_raw = gdat_data.error_bar;
gdat_data.error_bar = gdat_data.error_bar_raw * diag(tracefirrat_data);
gdat_data.data_fullpath=[gdat_data.data_fullpath '; fir_thom_rat ; _raw without 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