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
%
% Note for liuqe parameter: At this time we have liuqe_fortran (which used to be the default) and liuqe_matlab (which is new and now becomes/is the default)
% We still have liuqe1, 2 and 3 nodes for either of them, so you can choose:
% 'liuqe',1 (2 or 3): to get the default liuqe 1, 2 or 3 (which is now the matlab version)
% 'liuqe',11 (12 or 13): to get liuqe_fortran nodes 1, 2 or 3
% 'liuqe',21 (22 or 23): to get liuqe_matlab nodes 1, 2 or 3 (may be needed if we set the default to liuqe_fortran for old shots)
% 'liuqe',-1 : to get FBTE result
%
% 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)
% a7 = gdat(48836,'static("r_m")[$1]'',''001'); % note first and last quote of tdi argument added by gdat
% a7 = gdat(48836,'tcv_eq(''''psi'''',''''liuqe.m'''')'); % do not use double quote char so 'liuqe',11 will work
% a8 = gdat(64770,['\magnetics::bpol_003[*,$1]'',{''001''; ''030''}']); one string full expression
% a9 = gdat(64770,{'\magnetics::vloop[*,$1]',{'''001''','''A_001'''}}); cell array (other option)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Remote data access for TCV
% You need to make a "tunnel" in one unix session from the remote node using for example:
% ssh -L 5555:tcvdata.epfl.ch:8000 sauter@lac911.epfl.ch % to create a tunnel to port 5555
%
% Then in another unix session on the remote node, in matlab and with the appropriate mds path do:
% >> mdsconnect('localhost:5555')
% >> mdsvalue('1+2') % should return 3 if correctly connected
%
%
% 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;
gdat_params.trialindx = [];
% 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; % do not use xlsx anymore but scan respective machine_requests_mapping.m files
data_request_names = get_data_request_names_from_gdat_xxx(default_machine);
% 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 = sort(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_eff = data_request;
data_request = data_request; % should not lower until see if keyword
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 or not shot asked for
try
if isempty(mdsipmex(8))
shot_mds = shot;
else
shot_mds = mdsipmex(2,'$SHOT');
end
catch
shot_mds = shot;
end
if isnumeric(shot_mds); shot = shot_mds; end
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')
if ~strcmp(data_request.data_request,'ids')
warning(['shot cannot be opened with ' data_request.data_request]);
return
end
elseif ischar(data_request)
if ~strcmp(data_request,'ids')
warning(['shot cannot be opened with ' data_request]);
return
end
else
warning(['shot cannot be opened']);
return
end
end
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_eff = data_request.data_request;
% some request are case sensitive because mds names can be like \magnetics::ipol[*,"OH_001"]
data_request.data_request = 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:
i_liuqe_set_in_pairs = 0; % to know if liuqe was not set explicitely thus use value in expression later
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})
if strcmp(varargin_eff{i},'liuqe')
i_liuqe_set_in_pairs = 1;
end
% enforce lower case for any character driven input
if ischar(varargin_eff{i+1}) && ~strcmp('path',varargin_eff{i})
gdat_params.(lower(varargin_eff{i})) = varargin_eff{i+1}; % cannot enforce lower(params value) if filename or else are case sensitive
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
% if it is a request_keyword copy it:
if ischar(data_request_eff) || length(data_request_eff)==1
ij=strmatch(lower(data_request_eff),data_request_names_all,'exact');
else
ij=[];
end
if ~isempty(ij);
gdat_data.gdat_request = data_request_names_all{ij};
gdat_params.data_request = gdat_data.gdat_request;
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; % This means that before here it is gdat_params which should be updated
% 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
liuqe_matlab = 1; % now default should be matlab liuqe nodes
if liuqe_version<0 || (liuqe_version > 10 && liuqe_version < 20)
liuqe_matlab = 0;
end
liuqe_version_eff = mod(liuqe_version,10);
substr_liuqe = '';
substr_liuqe_tcv_eq = '';
if liuqe_version_eff==2 || liuqe_version_eff==3 || (liuqe_matlab==1 && liuqe_version_eff==1)
substr_liuqe = ['_' num2str(liuqe_version_eff)];
end
if liuqe_version_eff==2 || liuqe_version_eff==3
substr_liuqe_tcv_eq = num2str(liuqe_version_eff);
end
mapping_for_tcv = tcv_requests_mapping(data_request_eff,shot);
gdat_data.label = mapping_for_tcv.label;
% special treatment for model shot=-1 or preparation shot >=100'000
begstr = '';
if (iscell(mapping_for_tcv.expression) || isempty(strfind(mapping_for_tcv.expression,'\rtc::'))) && ...
~isempty(shot) && (shot==-1 || (shot>=100000 && shot < 200000) || liuqe_version==-1 )
% requires FBTE
liuqe_version_eff = -1;
begstr = 'tcv_eq( "';
substr_liuqe = '", "FBTE" )';
end
% should replace all above by just psitbx_str...
liuqe_matlab = 1;
switch liuqe_version
case {-1}, liuqe_ext=''; psitbx_str='FBTE'; liuqe_matlab = 0;
case {1,21}, liuqe_ext=''; psitbx_str='LIUQE.M';
case {11}, liuqe_ext=''; psitbx_str='LIUQE';liuqe_matlab = 0;
case {2, 3, 22, 23}, liuqe_ext=['_' num2str(mod(liuqe_version,10))]; psitbx_str=['LIUQE.M' num2str(mod(liuqe_version,10))];
case {12,13}, liuqe_ext=['_' num2str(mod(liuqe_version,10))]; psitbx_str=['LIUQE' num2str(mod(liuqe_version,10))];liuqe_matlab = 0;
otherwise, error(['Unknown LIUQE version, liuqe = ' num2str(liuqe_version)]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specifications on how to get the data provided in tcv_requests_mapping
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_eff==-1
if ~iscell(mapping_for_tcv.expression) && any(strfind(mapping_for_tcv.expression,'\rtc::'))
mdsdisconnect;
ishot = mdsopen('rtc',shot); % sub-tree
if any(strfind(data_request_eff,'\rtc::'))
data_request_eff = data_request_eff(7:end);
end
mapping_for_tcv.expression = mapping_for_tcv.expression(7:end);
elseif shot==-1 || liuqe_version_eff==-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 liuqe_matlab==0 && strcmp(mapping_for_tcv.method,'tdiliuqe'); substr_tdi = substr_liuqe; end
if iscell(mapping_for_tcv.expression)
if length(mapping_for_tcv.expression)>0
aaa_tmp = regexp(mapping_for_tcv.expression{1},{'toray','conf','ibs','proffit','astra'});
if ~isempty(gdat_data.gdat_params.trialindx) && gdat_data.gdat_params.trialindx >= 0 ...
&& ~isempty([aaa_tmp{:}])
% contains(mapping_for_tcv.expression{1},{'toray','conf','ibs','proffit','astra'}) % OS: contains not available in matlab850
ij = findstr(mapping_for_tcv.expression{1},':trial');
if isempty(ij)
aa = [mapping_for_tcv.expression{1} ':trial'];
if mdsdata(['node_exists("' aa '")']) > 0
mapping_for_tcv.expression_orig = mapping_for_tcv.expression;
mapping_for_tcv.expression{1} = aa;
end
end
end
% series of arguments for tdi given in each cell
if liuqe_matlab==1
ij = findstr(mapping_for_tcv.expression{1},'equil_');
if ~isempty(ij)
mapping_for_tcv.expression{1}(ij+5:ij+6) = substr_liuqe;
end
end
eval_expr = ['tdi(''' mapping_for_tcv.expression{1} ''''];
for i=2:length(mapping_for_tcv.expression)
if ~iscell(mapping_for_tcv.expression{i}) || length(mapping_for_tcv.expression{i}) == 1
eval_expr = [eval_expr ',''' mapping_for_tcv.expression{i} ''''];
else
eval_expr = [eval_expr ',{' mapping_for_tcv.expression{i}{1}];
for j=2:length(mapping_for_tcv.expression{i})
eval_expr=[eval_expr ',' mapping_for_tcv.expression{i}{j}];
end
eval_expr = [eval_expr '}'];
end
end
eval_expr = [eval_expr ');'];
aatmp = eval(eval_expr);
else
% empty or wrong expression
error_status=701;
return
end
else
do_liuqem2liuqef = 0;
if liuqe_version_eff==-1
% trialindx not relevant for fbte
mapping_for_tcv_expression_eff = mapping_for_tcv.expression;
if length(mapping_for_tcv.expression)>8 && strcmp(lower(mapping_for_tcv.expression(1:8)),'\results')
mapping_for_tcv_expression_eff = mapping_for_tcv.expression(11:end);
elseif findstr('tcv_eq',lower(mapping_for_tcv.expression))
ij = regexpi(mapping_for_tcv.expression,'''''');
if length(ij)<2
disp(['expected more double quotes in mapping_for_tcv.expression = ' mapping_for_tcv.expression]);
return
else
mapping_for_tcv_expression_eff = mapping_for_tcv.expression(ij(1)+2:ij(2)-1);
end
end
eval_expr = ['tdi(''' begstr mapping_for_tcv_expression_eff substr_liuqe ''');']
else
aaa_tmp = regexp(mapping_for_tcv.expression,{'toray','conf','ibs','proffit','astra'});
if ~isempty(gdat_data.gdat_params.trialindx) && gdat_data.gdat_params.trialindx >= 0 ...
&& ~isempty([aaa_tmp{:}])
% contains(mapping_for_tcv.expression,{'toray','conf','ibs','proffit','astra'}) % OS: not available in matlab850
ij = findstr(mapping_for_tcv.expression,':trial');
if isempty(ij)
aa = [mapping_for_tcv.expression ':trial'];
if mdsdata(['node_exists("' aa '")']) > 0 || mdsdata(['node_exists("\' aa '")']) > 0
mapping_for_tcv.expression_orig = mapping_for_tcv.expression;
mapping_for_tcv.expression = aa;
end
end
end
if liuqe_matlab==1
ij = findstr(mapping_for_tcv.expression,'equil_');
if ~isempty(ij)
mapping_for_tcv.expression(ij+5:ij+6) = substr_liuqe;
end
ij = regexpi(mapping_for_tcv.expression,'LIUQE\.M.?','once');
if ~isempty(ij)
ichar_after_liuqe = 7;
if strcmp(mapping_for_tcv.expression(ij+ichar_after_liuqe),'2') || ...
strcmp(mapping_for_tcv.expression(ij+ichar_after_liuqe),'3')
if i_liuqe_set_in_pairs==0
gdat_params.liuqe = str2num(mapping_for_tcv.expression(ij+ichar_after_liuqe));
gdat_data.gdat_params = gdat_params;
substr_liuqe_tcv_eq = mapping_for_tcv.expression(ij+ichar_after_liuqe);
end
ichar_after_liuqe = 8;
end
mapping_for_tcv.expression = [mapping_for_tcv.expression(1:ij+6) substr_liuqe_tcv_eq mapping_for_tcv.expression(ij+ichar_after_liuqe:end)];
else
% check if liuqe, liuqe2 or liuqe3 is given in expression
ij = regexpi(mapping_for_tcv.expression,'LIUQE[^\.]','once');
if ~isempty(ij)
ichar_after_liuqe = 5;
if strcmp(mapping_for_tcv.expression(ij+ichar_after_liuqe),'2') || ...
strcmp(mapping_for_tcv.expression(ij+ichar_after_liuqe),'3')
if i_liuqe_set_in_pairs==0
gdat_params.liuqe = 10+str2num(mapping_for_tcv.expression(ij+ichar_after_liuqe));
gdat_data.gdat_params = gdat_params;
substr_liuqe_tcv_eq = mapping_for_tcv.expression(ij+ichar_after_liuqe);
end
ichar_after_liuqe = 6;
else
if i_liuqe_set_in_pairs==0
gdat_params.liuqe = 11;
gdat_data.gdat_params = gdat_params;
substr_liuqe_tcv_eq = '';
end
end
if i_liuqe_set_in_pairs==1 && liuqe_matlab==1
% enforce matlab liuqe version even matlab if asked for
substr_liuqe_tcv_eq = ['.M' substr_liuqe_tcv_eq];
end
mapping_for_tcv.expression = [mapping_for_tcv.expression(1:ij+4) substr_liuqe_tcv_eq mapping_for_tcv.expression(ij+ichar_after_liuqe:end)];
end
end
else
ij = regexpi(mapping_for_tcv.expression,'LIUQE\.M.?','once');
if ~isempty(ij)
mapping_for_tcv.expression = [mapping_for_tcv.expression(1:ij+4) substr_liuqe_tcv_eq mapping_for_tcv.expression(ij+7:end)];
do_liuqem2liuqef = 1;
end
end
eval_expr = ['tdi(''' mapping_for_tcv.expression ''');'];
end
% special cases for matching liuqe fortran and matlab
if liuqe_matlab == 0 && do_liuqem2liuqef==1
% assume tcv_eq(''keyword'',''LIUQE..'')
liuqe_fortran_matlab = liuqefortran2liuqematlab;
ij = regexpi(eval_expr,'''''');
if numel(ij)>=2
liuqe_keyword = eval_expr(ij(1)+2:ij(2)-1);
end
ij_row=strmatch(liuqe_keyword,liuqe_fortran_matlab(:,2),'exact');
if ~isempty(ij_row)
eval_expr = [eval_expr(1:ij(1)+1) liuqe_fortran_matlab{ij_row,1} eval_expr(ij(1)+2+length(liuqe_keyword):end)];
end
end
aatmp=eval(eval_expr);
if ~isempty(gdat_data.gdat_params.trialindx) && gdat_data.gdat_params.trialindx >= 0
if isfield(mapping_for_tcv,'expression_orig')
mapping_for_tcv.expression = mapping_for_tcv.expression_orig;
mapping_for_tcv = rmfield(mapping_for_tcv,'expression_orig');
% reduce data to only desired trialindx
nbdims = length(aatmp.dim);
bbtmp = aatmp;
aatmp.dim = bbtmp.dim(1:nbdims-1);
aatmp.dimunits = bbtmp.dimunits(1:nbdims-1);
abc(1:nbdims-1)={':,'};
strdata=sprintf('%s',abc{:});
eval(['aatmp.data = bbtmp.data(' strdata num2str(gdat_data.gdat_params.trialindx+1) ');']);
aatmp.help = [bbtmp.help ' :trial(...,trialindx=' num2str(gdat_data.gdat_params.trialindx) ')'];
end
end
end
if isempty(aatmp.data) || (isempty(aatmp.dim) && ischar(aatmp.data) && ~isempty(strfind(lower(aatmp.data),'no data')))% || 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 (nbdims>1 && size(gdat_data.data,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
if length(gdat_data.dim)>=mapping_for_tcv.timedim && mapping_for_tcv.timedim>0;
gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
mapping_for_tcv.gdat_timedim = mapping_for_tcv.timedim;
gdat_data.mapping_for.tcv = mapping_for_tcv;
end
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
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];
if iscell(mapping_for_tcv.expression)
gdat_data.label=[mapping_for_tcv.expression{1}];
else
gdat_data.label=[mapping_for_tcv.expression];
end
gdat_data.help = aatmp.help;
% 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 = setdiff(fieldnames(gdat_tmp),{'gdat_request','label'}); % could/should also remove label in any case
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;
if isfield(gdat_tmp,'help')
gdat_data.help = gdat_tmp.help;
else
gdat_data.help = [];
end
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
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','r_geom','a_minor_rho','r_geom_rho','rgeom_rho'}
if strcmp(data_request_eff,'r_geom'); data_request_eff = 'rgeom'; end
if strcmp(data_request_eff,'r_geom_rho'); data_request_eff = 'rgeom_rho'; end
% compute average minor or major radius (on z=zaxis normally)
if liuqe_matlab==0
nodenameeff=['tcv_eq(''r_max_psi'',''LIUQE' substr_liuqe_tcv_eq ''')'];
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=['tcv_eq(''r_min_psi'',''LIUQE' substr_liuqe_tcv_eq ''')'];
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); % points outside TCV vessel
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,'a_minor_rho')
gdat_data.data=0.5.*(rmaxpsi.data(:,:) - rminpsi.data(:,:));
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'];
elseif strcmp(data_request_eff,'rgeom_rho')
gdat_data.data=0.5.*(rmaxpsi.data(:,:) + rminpsi.data(:,:));
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
if strcmp(data_request_eff(end-3:end),'_rho')
gdat_data.dim = rmaxpsi.dim;
gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
gdat_data.x = gdat_data.dim{1};
gdat_data.dimunits = rmaxpsi.dimunits(2);
else
gdat_data.dim = rmaxpsi.dim(2);
gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
gdat_data.dimunits = rmaxpsi.dimunits(2);
end
if any(strcmp(fieldnames(rmaxpsi),'units'))
gdat_data.units = rmaxpsi.units;
end
else
if isempty(substr_liuqe); substr_liuqe = '_1'; end
if strcmp(data_request_eff,'a_minor') % to update when ready, still buggy, this should be rho,t
nodenameeff=['tcv_eq(''a_minor_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
elseif strcmp(data_request_eff,'a_minor_rho')
nodenameeff=['tcv_eq(''r_out_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
nodenameeff2=['tcv_eq(''r_in_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
nodenameeff=['tcv_eq(''r_out'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
nodenameeff2=['tcv_eq(''r_in'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
elseif strcmp(data_request_eff,'rgeom')
nodenameeff=['tcv_eq(''r_geom_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
elseif strcmp(data_request_eff,'rgeom_rho')
nodenameeff=['tcv_eq(''r_out_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
nodenameeff2=['tcv_eq(''r_in_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
nodenameeff=['tcv_eq(''r_out'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
nodenameeff2=['tcv_eq(''r_in'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
else
if (gdat_params.nverbose>=1); warning(['should not be in this case with data_request_eff = ' data_request_eff]); end
return
end
aatmp = tdi(nodenameeff);
if strcmp(data_request_eff(end-3:end),'_rho')
aatmp2 = tdi(nodenameeff2);
if strcmp(data_request_eff,'a_minor_rho')
gdat_data.data = 0.5.*(aatmp.data-aatmp2.data);
elseif strcmp(data_request_eff,'rgeom_rho')
gdat_data.data = 0.5.*(aatmp.data+aatmp2.data);
end
gdat_data.dim = aatmp.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes
gdat_data.x = gdat_data.dim{1};
gdat_data.data_fullpath=[nodenameeff ' and ' nodenameeff2];
else
gdat_data.data = aatmp.data(end,:)';
gdat_data.dim = aatmp.dim(2);
aatmp.dimunits = aatmp.dimunits(2);
gdat_data.data_fullpath=[nodenameeff];
end
gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
gdat_data.units = aatmp.units(end);
gdat_data.dimunits = aatmp.dimunits;
end
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'zgeom','z_geom'}
% compute average minor or major radius (on z=zaxis normally)
%
data_request_eff = 'zgeom';
gdat_data.gdat_request = data_request_eff;
gdat_data.gdat_params.gdat_request = gdat_data.gdat_request;
gdat_params = gdat_data.gdat_params;
if liuqe_matlab==0
nodenameeff=['tcv_eq(''z_contour'',''LIUQE' substr_liuqe_tcv_eq ''')'];
else
if isempty(substr_liuqe); substr_liuqe = '_1'; end
nodenameeff=['tcv_eq(''z_edge'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
end
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
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'b0'}
% B0 at R0=0.88
r0exp=0.88;
if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source) ...
&& length(gdat_data.gdat_params.source)>=5 && strcmp(lower(gdat_data.gdat_params.source(1:5)),'liuqe')
% expect liuqe or liuqe.m to have liuqe time-base, otherwise give full time base
else
gdat_data.gdat_params.source = 'iphi';
end
if liuqe_version_eff==-1
nodenameeff = 'tcv_eq("BZERO","FBTE")';
tracetdi=tdi(nodenameeff);
gdat_data.data = tracetdi.data;
else
if strcmp(lower(gdat_data.gdat_params.source),'iphi')
nodenameeff=['\magnetics::iphi'];
tracetdi=tdi(nodenameeff);
added_correction = 0.996; % correction already in liuqe.f
added_correction_str = ['96*mu0/2piR0 * ' num2str(added_correction) ' * '];
gdat_data.data=192.E-07 * added_correction *tracetdi.data/r0exp;
else
if liuqe_matlab==0
added_correction = 1.0; % correction already in liuqe.f
added_correction_str = [''];
nodenameeff = ['tcv_eq(''BZERO'',''LIUQE' substr_liuqe_tcv_eq ''')'];
else
if isempty(substr_liuqe); substr_liuqe = '_1'; end
nodenameeff=['tcv_eq(''BZERO'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
added_correction = 0.996; % correction removed in liuqe.m at this stage
added_correction_str = [num2str(added_correction) ' * '];
end
tracetdi=tdi(nodenameeff);
gdat_data.data = tracetdi.data .* added_correction;
end
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=[added_correction_str 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'];
gdat_data.r0 = r0exp;
% At this stage implement time_out at resulting structure as post processing, to allow various interpolation and other options.
% Should implement at server level at some point, could exend tdi_os
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'betan', 'beta_tor_norm'}
% 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;
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
% should already have all data at times time_out since called gdat_tcv above
end
ij=find(isfinite(ip.data));
% use interpos linear with option 21 to have cst extrapolation
if numel(ij) > 1
ip_t = interpos(21,ip.dim{1}(ij),ip.data(ij),gdat_data.t);
elseif numel(ij) == 1
ip_t = ip.data(ij);
else
warning('no valid data for ip, cannot fit')
return
end
ij=find(isfinite(b0.data));
if numel(ij) > 1
b0_t = interpos(21,b0.dim{1}(ij),b0.data(ij),gdat_data.t);
elseif numel(ij) == 1
b0_t = b0.data(ij);
else
warning('no valid data for b0, cannot fit')
return
end
ij=find(isfinite(a_minor.data));
if numel(ij) > 1
a_minor_t = interpos(21,a_minor.dim{1}(ij),a_minor.data(ij),gdat_data.t);
elseif numel(ij) == 1
a_minor_t = a_minor.data(ij);
else
warning('no valid data for a_minor, cannot fit')
return
end
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','cxrs_rho'}
%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;
gdat_data.cxrs_params.defaults = cxrs_params;
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
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,'cxrs_time_interval') && ~isempty(params_eff.cxrs_time_interval) && length(params_eff.cxrs_time_interval)>=2
cxrs_time_interval = params_eff.cxrs_time_interval;
if ~isfield(params_eff,'cxrs_plot')
cxrs_plot=1;
end
else
cxrs_time_interval = [];
end
gdat_data.gdat_params.cxrs_time_interval = cxrs_time_interval;
gdat_data.gdat_params.cxrs_plot = cxrs_plot;
if isfield(params_eff,'cxrs_xout')
cxrs_xout = params_eff.cxrs_xout;
else
cxrs_xout = [linspace(0,1.,201) [1.01:0.01:1.05]];
end
gdat_data.gdat_params.cxrs_xout = cxrs_xout;
cxrs_params.prof.all.xout = cxrs_xout;
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,cxrs_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.cxrs_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.cxrs_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.cxrs_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(cxrs_time_interval)
gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[' num2str(source) '],[],cxrs_params); % with cxrs_params'];
else
gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[' num2str(source) '],[' num2str(cxrs_time_interval) '],cxrs_params); % with cxrs_params'];
end
gdat_data.dimunits{1} = '';
gdat_data.dimunits{2} = 's';
% add grids_1d to have rhotor, etc if cxrs_rho was asked for
if strcmp(data_request_eff,'cxrs_rho')
gdat_data = get_grids_1d(gdat_data,2,1,gdat_params.nverbose);
else
gdat_data = get_grids_1d(gdat_data,2,0,gdat_params.nverbose);
end
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) && ~isempty(gdat_data.t) ...
&& ~isempty(gdat_data.data) && isnumeric(gdat_data.data)
others_as_data={'error_bar'};
size_data = size(gdat_data.data);
for i=1:length(sub_nodes_out)
fn_sub = fieldnames(gdat_data.(sub_nodes_out{i}));
fn_sub = setdiff(fn_sub,'units');
for j=1:length(fn_sub)
fn_sub_sub=fieldnames(gdat_data.(sub_nodes_out{i}).(fn_sub{j}));
for k=1:length(fn_sub_sub)
if ~isempty(gdat_data.(sub_nodes_out{i}).(fn_sub{j}).(fn_sub_sub{k}))
try
size_leaf = size(gdat_data.(sub_nodes_out{i}).(fn_sub{j}).(fn_sub_sub{k}));
if size_data(gdat_data.mapping_for.tcv.timedim) == size_leaf(gdat_data.mapping_for.tcv.timedim) ...
&& length(size_data)==length(size_leaf)
others_as_data{end+1} = [sub_nodes_out{i} '.' fn_sub{j} '.' fn_sub_sub{k}];
end
catch
end
end
end
end
end
if gdat_data.gdat_params.nverbose >= 3
disp(['in cxrs, fields with data reduced to match time_out:'])
others_as_data
end
[gdat_data] = gdat2time_out(gdat_data,1,others_as_data);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'eqdsk'}
%
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time)
disp('time_out parameter not used for eqdsk, uses time')
gdat_data.gdat_params = rmfield(gdat_data.gdat_params,'time_out');
else
warning('time_out parameter not relevant for eqdsk, use time to give time array for outputs')
return;
end
end
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
default_write_eqdsk = 1;
if isfield(gdat_data.gdat_params,'write') && ~isempty(gdat_data.gdat_params.write)
if ischar(gdat_data.gdat_params.write)
if strcmp(lower(gdat_data.gdat_params.write),'no')
gdat_data.gdat_params.write = 0;
elseif strcmp(lower(gdat_data.gdat_params.write),'yes')
gdat_data.gdat_params.write = 1;
else
if (gdat_params.nverbose>=3); disp(['expects 0 or 1, event yes or no for write option, not: ' ...
gdat_data.gdat_params.write ', use default 1']);
end
gdat_data.gdat_params.write = [];
end
end
if gdat_data.gdat_params.write~=1 && gdat_data.gdat_params.write~=0
gdat_data.gdat_params.write = default_write_eqdsk;
end
else
gdat_data.gdat_params.write = default_write_eqdsk;
end
default_map_eqdsk_psirz = 0;
if isfield(gdat_data.gdat_params,'map_eqdsk_psirz') && ~isempty(gdat_data.gdat_params.map_eqdsk_psirz)
if ~isnumeric(gdat_data.gdat_params.map_eqdsk_psirz) || (gdat_data.gdat_params.map_eqdsk_psirz~=1 && gdat_data.gdat_params.map_eqdsk_psirz~=0)
gdat_data.gdat_params.map_eqdsk_psirz = default_map_eqdsk_psirz_eqdsk;
end
else
gdat_data.gdat_params.map_eqdsk_psirz = default_map_eqdsk_psirz;
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;
nrz_out = [129,129];
if isfield(gdat_data.gdat_params,'nrz_out') && ~isempty(gdat_data.gdat_params.nrz_out)
nrz_out = gdat_data.gdat_params.nrz_out;
else
gdat_data.gdat_params.nrz_out = nrz_out;
end
gdat_data.gdat_params.nrz_out = nrz_out;
sources_available = {'liuqe','chease'};
if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source) && ~isempty(intersect(sources_available,lower(gdat_data.gdat_params.source)))
source = gdat_data.gdat_params.source;
else
gdat_data.gdat_params.source = 'liuqe';
end
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,nrz_out);
cocos_read_results_for_chease = 17;
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
ii = regexpi(fnames_readresults{4},'eqdsksigns');
if ~isempty(ii)
eqdskval=read_eqdsk(fnames_readresults{4},cocos_read_results_for_chease,0,[],[],1); % read_results now has relevant cocos LIUQE= 17
else
disp('wrong name check order of eqdsk outputs')
return
end
for i=1:length(fnames_readresults)
unix(['rm ' fnames_readresults{i}]);
end
%
% run CHEASE if asked
%
cocos_in = 2;
if strcmp(lower(gdat_data.gdat_params.source),'chease')
% will give cocos_out (becoming cocos_in)=2 by default
[fname_out,globalsvalues,namelist_struct,namelistfile_eff] = run_chease(2,eqdskval,cocos_read_results_for_chease);
ij = strfind(fname_out,'EQDSK_COCOS_02.OUT');
i = [];
for i=1:length(ij)
if ~isempty(ij(i)); break;end
end
eqdsk_cocos_in = read_eqdsk(fname_out{i},2,0);
else
%
% transform to cocos=2 since read_results originally assumed it was cocos=2
[eqdsk_cocos_in, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdskval,[cocos_read_results_for_chease cocos_in]);
end
%
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)
if gdat_data.gdat_params.write==1
write_eqdsk(fnamefull,eqdsk_cocos_in,cocos_in,[],[],[],1);
end
% 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
if gdat_data.gdat_params.write==1
gdat_data.eqdsk(itime) = write_eqdsk(fnamefull,eqdsk_cocosout,cocos_out);
else
gdat_data.eqdsk(itime) = eqdsk_cocosout;
end
if gdat_data.gdat_params.map_eqdsk_psirz==1
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
% do not map all psi(r,z) onto same mesh and leave .data, .dim empty (much faster)
end
else
if gdat_data.gdat_params.write==1
gdat_data.eqdsk = write_eqdsk(fnamefull,eqdsk_cocosout,cocos_out);
else
gdat_data.eqdsk = eqdsk_cocosout;
end
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);
if gdat_data.gdat_params.map_eqdsk_psirz==1
gdat_data.data_fullpath=['psi(R,Z,t) on same R,Zmesh in .data and eqdsk(itime) from read_eqdsk from LIUQE' num2str(liuqe_version) ';zshift=' num2str(zshift)];
else
gdat_data.data_fullpath=['eqdsk(itime) from read_eqdsk from LIUQE' num2str(liuqe_version) ';zshift=' num2str(zshift)];
end
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 {'gas', 'gas_flux', 'gas_request', 'gas_feedforward','gas_valve'}
params_eff = gdat_data.gdat_params;
params_eff.data_request = '\diagz::flux_gaz:piezo_1:flux';
gasflux = gdat_tcv(gdat_data.shot,params_eff);
gdat_data.gas_flux = gasflux;
params_eff.data_request = '\hybrid::mat_m_signals:output_019';
gasrequest = gdat_tcv(gdat_data.shot,params_eff);
gdat_data.gas_request_volt = gasrequest;
%volt_to_mlpers = max(0.,(in_volt-0.6879)*(72.41/(4.2673-0.6879))
gdat_data.gas_request_flux = gasrequest;
gdat_data.gas_request_flux.data = max(0.,72.41.*(gasrequest.data-0.6879)./(4.2673-0.6879));
gdat_data.gas_request_flux.units = gasflux.units;
params_eff.data_request = '\draw_feedfor_gas:alim_001';
gasrequest_ref = gdat_tcv(gdat_data.shot,params_eff); gasrequest_ref.units = 'V';
gdat_data.gas_feedforward_volt = gasrequest_ref;
gdat_data.gas_feedforward_flux = gasrequest_ref;
gdat_data.gas_feedforward_flux.data = max(0.,72.41.*(gasrequest_ref.data-0.6879)./(4.2673-0.6879));
gdat_data.label = data_request_eff;
switch data_request_eff
case {'gas', 'gas_flux'}
gdat_data.data = gdat_data.gas_flux.data;
gdat_data.units = gdat_data.gas_flux.units;
gdat_data.t = gdat_data.gas_flux.t;
gdat_data.data_fullpath = gdat_data.gas_flux.data_fullpath;
case {'gas_request'}
gdat_data.data = gdat_data.gas_request_volt.data;
gdat_data.units = gdat_data.gas_request_volt.units;
gdat_data.t = gdat_data.gas_request_volt.t;
gdat_data.data_fullpath = gdat_data.gas_request_volt.data_fullpath;
case {'gas_feedforward'}
gdat_data.data = gdat_data.gas_feedforward_volt.data;
gdat_data.units = gdat_data.gas_feedforward_volt.units;
gdat_data.t = gdat_data.gas_feedforward_volt.t;
gdat_data.data_fullpath = gdat_data.gas_feedforward_volt.data_fullpath;
case {'gas_valve'}
if ~isfield(params_eff,'source') || isempty(params_eff.source)
params_eff.source = {{1,'D2'}, {3,'N2'}};
end
if iscell(params_eff.source) && numel(params_eff.source)>0 && numel(params_eff.source{1})>1
gdat_data.units = '';
for i=1:numel(params_eff.source)
try
[aat,aad] = get_gas_flux(shot,params_eff.source{i}{1},params_eff.source{i}{2});
gdat_data.gas_valve(i).t = aat;
gdat_data.gas_valve(i).data = aad;
gdat_data.gas_valve(i).valve_nb = params_eff.source{i}{1};
gdat_data.gas_valve(i).gas_type = params_eff.source{i}{2};
gdat_data.gas_valve(i).units = 'Molecules/s';
if ~strcmp(gdat_data.gas_valve(i).gas_type(end),'2'); gdat_data.gas_valve(i).units = 'atoms/s'; end
gdat_data.data(:,i) = gdat_data.gas_valve(i).data;
gdat_data.t = gdat_data.gas_valve(i).t;
gdat_data.dim{2}(i) = gdat_data.gas_valve(i).valve_nb;
gdat_data.units{i} = gdat_data.gas_valve(i).units;
gdat_data.label = [gdat_data.label ',' gdat_data.gas_valve(i).gas_type];
if i==1
gdat_data.dim{1} = gdat_data.t;
gdat_data.dimunits{1} = 's';
gdat_data.dimunits{2} = 'valve nb';
gdat_data.mapping_for.tcv.gdat_timedim = 1;
gdat_data.mapping_for.tcv.timedim = 1;
gdat_data.data_fullpath = 'using get_gas_flux(shot,''gdat_params.source{1}'',''gdat_params.source{2}'')';
end
catch ME
gdat_data.gas_valve(i).t = [];
gdat_data.gas_valve(i).data = [];
gdat_data.gas_valve(i).valve_nb = params_eff.source{i}{1};
gdat_data.gas_valve(i).gas_type = params_eff.source{i}{2};
gdat_data.gas_valve(i).units = '';
end
end
end
otherwise
error('gas option not defined')
end
if ~strcmp(data_request_eff,'gas_valve')
gdat_data.dim{1} = gdat_data.t;
gdat_data.dimunits{1} = 's';
end
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'halphas'}
channels = [1:18];
if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels)
channels = gdat_data.gdat_params.channels;
end
% '\base::pd:pd_001';
params_eff = gdat_data.gdat_params;
for i=1:length(channels)
params_eff.data_request=['\base::pd:pd_' num2str(channels(i),'%.3d')];
try
halphas=gdat_tcv(gdat_data.shot,params_eff); % note: no need to set .doplot=0 since gdat_tcv does not call gdat_plot in any case
if i == 1;
gdat_data.data = NaN*ones(max(channels),length(halphas.t));
gdat_data.t = halphas.t;
end
gdat_data.data(channels(i),:) = halphas.data;
catch
disp([params_eff.data_request ' not ok']);
end
end
gdat_data.label = ['\base::pd:pd_' num2str(min(channels),'%.3d') ' to ' num2str(max(channels),'%.3d')];
gdat_data.x = channels;
gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits{1} = 'channel #';
gdat_data.dimunits{2} = 's';
gdat_data.mapping_for.tcv.gdat_timedim = 2;
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ids'}
ids_empty_path = fullfile(fileparts(mfilename('fullpath')),'..','TCV_IMAS','ids_empty');
params_eff = gdat_data.gdat_params;
if isfield(params_eff,'source') && ~isempty(params_eff.source)
ids_top_name = params_eff.source;
else
ids_top_name = [];
warning('gdat:EmptyIDSName','Need an ids name in ''source'' parameter\n check substructure gdat_params.sources_available for an ids list');
addpath(ids_empty_path);
assert(~~exist('ids_list_all','file'),'could not find ids_list_all.m in %s',ids_empty_path);
gdat_data.gdat_params.sources_available = ids_list_all; % avoid adding ids_list.m which would override ids_list from IMAS
rmpath(ids_empty_path);
return
end
ids_gen_ok = ~~exist('ids_gen','file');
if ids_gen_ok
ids_empty = ids_gen(ids_top_name); % generate ids from gateway function ids_gen
else
% load empty ids structure from template file
fname = sprintf('ids_empty_%s',ids_top_name);
try
assert(~~exist(ids_empty_path,'dir'),'folder %s not found',ids_empty_path);
addpath(ids_empty_path);
assert(~~exist(fname,'file'),'file %s does not exist in %s',fname,ids_empty_path);
ids_empty = eval(fname);
disp(['use structure in .mat file: ' ids_empty_path '/' fname]);
rmpath(ids_empty_path);
catch ME
fprintf('Could not load empty template for %s\n',ids_top_name);
fprintf('Available templates:\n');
disp(dir(ids_empty_path));
rmpath(ids_empty_path);
rethrow(ME);
end
end
if ~isfield(gdat_data.gdat_params,'error_bar') || isempty(gdat_data.gdat_params.error_bar)
gdat_data.gdat_params.error_bar = 'delta';
end
if ~isfield(gdat_data.gdat_params,'cocos_in') || isempty(gdat_data.gdat_params.cocos_in)
gdat_data.gdat_params.cocos_in = 17;
end
if ~isfield(gdat_data.gdat_params,'cocos_out') || isempty(gdat_data.gdat_params.cocos_out)
gdat_data.gdat_params.cocos_out = 11;
end
if ~isfield(gdat_data.gdat_params,'ipsign_out') || isempty(gdat_data.gdat_params.ipsign_out)
gdat_data.gdat_params.ipsign_out = 0;
end
if ~isfield(gdat_data.gdat_params,'b0sign_out') || isempty(gdat_data.gdat_params.b0sign_out)
gdat_data.gdat_params.b0sign_out = 0;
end
if ~isfield(gdat_data.gdat_params,'ipsign_in') || isempty(gdat_data.gdat_params.ipsign_in)
if gdat_data.gdat_params.ipsign_out~=0
gdat_data.gdat_params.ipsign_in = mdsvalue('\pcs::data:iohfb');
else
gdat_data.gdat_params.ipsign_in = 0;
end
end
if ~isfield(gdat_data.gdat_params,'b0sign_in') || isempty(gdat_data.gdat_params.b0sign_in)
if gdat_data.gdat_params.b0sign_out~=0
gdat_data.gdat_params.b0sign_in = mdsvalue('\pcs::data:if36fb');
else
gdat_data.gdat_params.b0sign_in = 0;
end
end
if strcmp(ids_top_name,'summary')
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
warning(['For summary ids, the present convention is to have time defined within get_ids_summary, but since time_out is provided,', ...
' it is not modified but there might be extra points']);
end
end
try
if ~isempty(shot)
[ids_top,ids_top_description] = feval(['tcv_get_ids_' ids_top_name],shot,ids_empty,gdat_data.gdat_params);
gdat_data.(ids_top_name) = ids_top;
gdat_data.([ids_top_name '_description']) = ids_top_description;
else
gdat_data.(ids_top_name) = ids_empty;
gdat_data.([ids_top_name '_description']) = ['shot empty so return default empty structure for ' ids_top_name];
end
catch ME_tcv_get_ids
disp(['there is a problem with: tcv_get_ids_' ids_top_name ...
' , may be check if it exists in your path or test it by itself'])
gdat_data.(ids_top_name) = ids_empty;
gdat_data.([ids_top_name '_description']) = getReport(ME_tcv_get_ids);
rethrow(ME_tcv_get_ids)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ec_data', 'aux', 'h_cd', 'nbi_data', 'ic_data', 'lh_data','ohm_data', 'bs_data'}
% note: same time array for all at main.data level, then individual at .ec, .nbi levels
% At this stage fill just eccd, later add nbi
sources_avail = {'ec','nbi','ohm','bs'}; % can be set in parameter source
% create empty structures for all, so in return one always have same substructres
for i=1:length(sources_avail)
gdat_data.(sources_avail{i}).data = [];
gdat_data.(sources_avail{i}).units = [];
gdat_data.(sources_avail{i}).dim=[];
gdat_data.(sources_avail{i}).dimunits=[];
gdat_data.(sources_avail{i}).t=[];
gdat_data.(sources_avail{i}).x=[];
gdat_data.(sources_avail{i}).data_fullpath=[];
gdat_data.(sources_avail{i}).label=[];
end
if ~isempty(gdat_data.gdat_params.trialindx) && gdat_data.gdat_params.trialindx < 0
gdat_data.gdat_params.trialindx = [];
end
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
switch data_request_eff
case 'ec_data'
gdat_data.gdat_params.source = {'ec'};
case 'ohm_data'
gdat_data.gdat_params.source = {'ohm'};
case 'bs_data'
gdat_data.gdat_params.source = {'bs'};
otherwise
gdat_data.gdat_params.source = sources_avail;
end
elseif ~iscell(gdat_data.gdat_params.source)
if ischar(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source);
if ~any(strmatch(gdat_data.gdat_params.source,lower(sources_avail)))
if (gdat_params.nverbose>=1)
warning(['source= ' gdat_data.gdat_params.source ' is not part of the available sources: ' sprintf('''%s'' ',sources_avail{:})]);
end
return
else
gdat_data.gdat_params.source = {gdat_data.gdat_params.source};
end
else
if (gdat_params.nverbose>=1); warning([' source parameter not compatible with: ' sprintf('''%s'' ',sources_avail{:})]); end
return
end
else
for i=1:length(gdat_data.gdat_params.source)
gdat_data.gdat_params.source{i} = lower(gdat_data.gdat_params.source{i});
if ~any(strmatch(gdat_data.gdat_params.source{i},lower(sources_avail)))
if gdat_data.gdat_params.nverbose>=1
warning(['source = ' gdat_data.gdat_params.source{i} ' not expected with data_request= ' data_request_eff])
end
end
end
end
% create structure for icd sources from params and complete with defaults
source_icd.ec = 'toray';
source_icd.nbi = '';
source_icd.ohm = 'ibs';
source_icd.bs = 'ibs';
for i=1:length(gdat_data.gdat_params.source)
if ~isfield(gdat_data.gdat_params,['source_' gdat_data.gdat_params.source{i}])
gdat_data.gdat_params.(['source_' gdat_data.gdat_params.source{i}]) = source_icd.(gdat_data.gdat_params.source{i});
else
source_icd.(gdat_data.gdat_params.source{i}) = gdat_data.gdat_params.(['source_' gdat_data.gdat_params.source{i}]);
end
end
mdsopen(shot);
field_for_main_data = 'cd_tot';
% add each source in main.data, on ohm time array
gdat_data.units = 'A';
gdat_data.label=[]; % label was defined in tcv_mapping_request as char so replace to allow cells
%
if any(strmatch('ec',gdat_data.gdat_params.source))
data_fullpath = '';
ec_help = '';
% EC
if strcmp(lower(source_icd.ec),'toray')
if isempty(gdat_data.gdat_params.trialindx)
[pabs_gyro,icdtot,pow_dens,currentdrive_dens,rho_dep_pow,drho_pow,pmax,icdmax,currentdrive_dens_w2,rho_dep_icd,drho_icd] = astra_tcv_EC_exp(shot); % centralized function for toray nodes
else
[pabs_gyro,icdtot,pow_dens,currentdrive_dens,rho_dep_pow,drho_pow,pmax,icdmax,currentdrive_dens_w2,rho_dep_icd,drho_icd] = astra_tcv_EC_exp(shot,[],[],[],[],[],gdat_data.gdat_params.trialindx); % centralized function for toray nodes
end
if gdat_data.mapping_for.tcv.gdat_timedim ==2
tgrid_to_change = {'pabs_gyro','icdtot','pow_dens','currentdrive_dens','rho_dep_pow','drho_pow','pmax', ...
'icdmax','currentdrive_dens_w2','rho_dep_icd','drho_icd'};
for i=1:length(tgrid_to_change)
eval([tgrid_to_change{i} '.tgrid = reshape(' tgrid_to_change{i} '.tgrid,1,numel(' tgrid_to_change{i} '.tgrid));']);
end
end
ec_help = 'from toray icdint with extracting of effective Icd for given launcher depending on nb rays used';
% All EC related quantities, each substructure should have at least fields data,x,t,units,dim,dimunits,label to be copied onto gdat_data
launchers_label = {'1','2','3','4','5','6','7','8','9','tot'};
launchers_grid = [1:10]';
% power deposition related:
ec_data.p_abs_plasma.data = pabs_gyro.data * 1e6;
ec_data.p_abs_plasma.data(end+1,:) = sum(ec_data.p_abs_plasma.data,1,'omitnan');
ec_data.p_abs_plasma.label = [strrep(pabs_gyro.comment,'MW','W') ' ; last index is total'];
ec_data.p_abs_plasma.units = 'W';
ec_data.p_abs_plasma.x = launchers_grid;
ec_data.p_abs_plasma.t =pabs_gyro.tgrid;
ec_data.p_abs_plasma.dim = {ec_data.p_abs_plasma.x, ec_data.p_abs_plasma.t};
ec_data.p_abs_plasma.dimunits = {launchers_label, 's'};
%
ec_data.p_dens.data = pow_dens.data * 1e6;
ec_data.p_dens.data(:,end+1,:) = sum(ec_data.p_dens.data,2,'omitnan');
ec_data.p_dens.label = [strrep(pow_dens.comment,'MW','W') ' ; last index is total'];
ec_data.p_dens.units = 'W/m^3';
ec_data.p_dens.x = pow_dens.rgrid';
ec_data.p_dens.rhotor_norm = ec_data.p_dens.x;
ec_data.p_dens.t = pow_dens.tgrid;
ec_data.p_dens.dim = {ec_data.p_dens.x, launchers_grid, ec_data.p_dens.t};
ec_data.p_dens.dimunits = {'rhotor_norm', launchers_label, 's'};
%
ec_data.max_pow_dens.data = pmax.data * 1e6;
ec_data.max_pow_dens.label = strrep(pmax.comment,'MW','W');
ec_data.max_pow_dens.units = 'W/m^3';
ec_data.max_pow_dens.x = [];
ec_data.max_pow_dens.t = pmax.tgrid;
ec_data.max_pow_dens.dim = {ec_data.max_pow_dens.t};
ec_data.max_pow_dens.dimunits = {'s'};
%
ec_data.rho_max_pow_dens.data = rho_dep_pow.data * 1e6;
ec_data.rho_max_pow_dens.label = strrep(rho_dep_pow.comment,'MW','W');
ec_data.rho_max_pow_dens.units = 'rhotor_norm';
ec_data.rho_max_pow_dens.x = [];
ec_data.rho_max_pow_dens.t = rho_dep_pow.tgrid;
ec_data.rho_max_pow_dens.dim = {ec_data.rho_max_pow_dens.t};
ec_data.rho_max_pow_dens.dimunits = {'s'};
%
ec_data.width_pow_dens.data = drho_pow.data;
ec_data.width_pow_dens.label = drho_pow.comment;
ec_data.width_pow_dens.units = 'rhotor_norm';
ec_data.width_pow_dens.x = [];
ec_data.width_pow_dens.t = drho_pow.tgrid;
ec_data.width_pow_dens.dim = {ec_data.width_pow_dens.t};
ec_data.width_pow_dens.dimunits = {'s'};
% current drive deposition related:
ec_data.cd_tot.data = icdtot.data * 1e6;
ec_data.cd_tot.data(end+1,:) = sum(ec_data.cd_tot.data,1,'omitnan');
ec_data.cd_tot.label = [strrep(icdtot.comment,'MA','A') ' ; last index is total'];
ec_data.cd_tot.units = 'A';
ec_data.cd_tot.x = launchers_grid;
ec_data.cd_tot.t = icdtot.tgrid;
ec_data.cd_tot.dim = {ec_data.cd_tot.x, ec_data.cd_tot.t};
ec_data.cd_tot.dimunits = {launchers_label, 's'};
%
ec_data.cd_dens.data = currentdrive_dens.data * 1e6;
ec_data.cd_dens.data(:,end+1,:) = sum(ec_data.cd_dens.data,2,'omitnan');
ec_data.cd_dens.label = [strrep(currentdrive_dens.comment,'MA','A') ' ; last index is total'];
ec_data.cd_dens.units = 'A/m^2';
ec_data.cd_dens.x = currentdrive_dens.rgrid';
ec_data.cd_dens.rhotor_norm = ec_data.cd_dens.x;
ec_data.cd_dens.t = currentdrive_dens.tgrid;
ec_data.cd_dens.dim = {ec_data.cd_dens.x, launchers_grid, ec_data.cd_dens.t};
ec_data.cd_dens.dimunits = {'rhotor_norm', launchers_label, 's'};
%
ec_data.max_cd_dens.data = icdmax.data * 1e6;
ec_data.max_cd_dens.label = strrep(icdmax.comment,'MA','A');
ec_data.max_cd_dens.units = 'A/m^2';
ec_data.max_cd_dens.x = [];
ec_data.max_cd_dens.t = icdmax.tgrid;
ec_data.max_cd_dens.dim = {ec_data.max_cd_dens.t};
ec_data.max_cd_dens.dimunits = {'s'};
%
ec_data.rho_max_cd_dens.data = rho_dep_icd.data;
ec_data.rho_max_cd_dens.label = rho_dep_icd.comment;
ec_data.rho_max_cd_dens.units = 'rhotor_norm';
ec_data.rho_max_cd_dens.x = [];
ec_data.rho_max_cd_dens.t = rho_dep_icd.tgrid;
ec_data.rho_max_cd_dens.dim = {ec_data.rho_max_cd_dens.t};
ec_data.rho_max_cd_dens.dimunits = {'s'};
%
ec_data.width_cd_dens.data = drho_icd.data;
ec_data.width_cd_dens.label = drho_icd.comment;
ec_data.width_cd_dens.units = 'rhotor_norm';
ec_data.width_cd_dens.x = [];
ec_data.width_cd_dens.t = drho_icd.tgrid;
ec_data.width_cd_dens.dim = {ec_data.width_cd_dens.t};
ec_data.width_cd_dens.dimunits = {'s'};
%
ec_data.cd_dens_doublewidth.data = currentdrive_dens_w2.data * 1e6;
ec_data.cd_dens_doublewidth.label = [strrep(currentdrive_dens_w2.comment,'MA','A') ' ; last index is total'];
for subfields={'x','rhotor_norm','t','dim','dimunits','units'}
ec_data.cd_dens_doublewidth.(subfields{1}) = ec_data.cd_dens.(subfields{1});
end
else
disp(['source_icd.ec = ' source_icd.ec ' not yet implemented, ask O. Sauter'])
ec_data.p_abs_plasma = [];
ec_data.p_abs_plasma(end+1,:) = [];
ec_data.p_abs_plasma_label = [];
ec_data.p_dens = [];
ec_data.p_dens(end+1,:) = [];
ec_data.p_dens_label = [];
ec_data.max_pow_dens = [];
ec_data.max_pow_dens_label = [];
ec_data.rho_max_pow_dens = [];
ec_data.rho_max_pow_dens_label = [];
ec_data.width_pow_dens = [];
ec_data.width_pow_dens_label = [];
% current drive deposition related:
ec_data.cd_tot = [];
ec_data.cd_tot(end+1,:) = [];
ec_data.cd_tot_label = [];
ec_data.cd_dens = [];
ec_data.cd_dens(:,end+1,:) = [];
ec_data.cd_dens_label = [];
ec_data.max_cd_dens = [];
ec_data.max_cd_dens_label = [];
ec_data.rho_max_cd_dens = [];
ec_data.rho_max_cd_dens_label = [];
ec_data.width_cd_dens = [];
ec_data.width_cd_dens_label = [];
ec_data.cd_dens_doublewidth = [];
ec_data.cd_dens_doublewidth_label = [];
ec_data.rho_tor_norm = [];
ec_data.t = [];
ec_data.launchers = [];
gdat_data.ec.ec_data = ec_data;
return
end
gdat_data.ec.ec_data = ec_data;
if isempty(icdtot.data) || isempty(icdtot.tgrid) || ischar(icdtot.data)
if (gdat_params.nverbose>=1)
warning(['problems loading data for ' source_icd.ec ...
' for data_request= ' data_request_eff]);
end
else
% now default is icdtot, will depend on request and data_out param of some kind
data_fullpath = ['from toray nodes using astra_tcv_EC_exp(shot), all results in .ec_data, subfield=' field_for_main_data ...
'in ec.data, .x, .t, .dim, .dimunits, .label, .units'];
for subfields={'data','x','t','units','dim','dimunits','label'}
gdat_data.ec.(subfields{1}) = gdat_data.ec.ec_data.(field_for_main_data).(subfields{1});
end
gdat_data.ec.data_fullpath = data_fullpath;
gdat_data.ec.help = ec_help;
% add to main, assume 1st one so just use this time base and x base
% should find launcher tot index
gdat_data.data(end+1,:) = gdat_data.ec.data(end,:);
gdat_data.t = gdat_data.ec.t;
if ischar(gdat_data.label); gdat_data.label = []; end; % label was defined in tcv_mapping_request as char so replace 1st time
gdat_data.label{end+1}=gdat_data.ec.label;
end
end
%
if any(strmatch('nb',gdat_data.gdat_params.source))
NBH_in_TCV = 0;
if shot >= 51641
% NBI
moderemote = mdsdata('data(\VSYSTEM::TCV_NBHDB_B["NBI:MODEREMOTE_FB"])');
lcs_mode = mdsdata('data(\VSYSTEM::TCV_NBHDB_I["NBI:LCS_MODE"])');
% NBH in TCV equiv moderemote='ON' AND lcs_mode = 9
NBH_in_TCV = strcmpi(strtrim(moderemote),'ON') && lcs_mode == 9;
else
% Nodes used in previous block only exist outside of Vista for shots after 51641
if any(shot == [51458 51459 51460 51461 51463 51465 51470 51472 ... % 29.JAN.2016
51628 51629 51631 51632 51633 ... % 09.FEB.2016
51639 51640 ... 51641 % 10.FEB.2016
])
NBH_in_TCV = 1;
end
end
if NBH_in_TCV
% should add reading from file at this stage ala summary Karpushov
end
end
%
if any(strmatch('ohm',gdat_data.gdat_params.source))
data_fullpath = '';
ohm_help = '';
if strcmp(lower(source_icd.ohm),'ibs')
ohm_data.cd_tot = gdat_tcv(shot,'\results::ibs:iohm');
if gdat_data.mapping_for.tcv.gdat_timedim ==2
tgrid_to_change = {'ohm_data.cd_tot.t','ohm_data.cd_tot.dim{1}'};
for i=1:length(tgrid_to_change)
eval([tgrid_to_change{i} ' = reshape(' tgrid_to_change{i} ',1,numel(' tgrid_to_change{i} '));']);
end
end
ohm_data.cd_tot.data = ohm_data.cd_tot.data';
ohm_data.cd_tot.units = strrep(ohm_data.cd_tot.units,'kA','A');
ohm_data.cd_dens = gdat_tcv(shot,'\results::ibs:johmav'); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R
if gdat_data.mapping_for.tcv.gdat_timedim ==2
tgrid_to_change = {'ohm_data.cd_dens.t','ohm_data.cd_dens.dim{2}'};
for i=1:length(tgrid_to_change)
eval([tgrid_to_change{i} ' = reshape(' tgrid_to_change{i} ',1,numel(' tgrid_to_change{i} '));']);
end
end
ohm_data.cd_dens.units = strrep(ohm_data.cd_dens.units,'kA','A');
abc=get_grids_1d(ohm_data.cd_dens,1,1);
ohm_data.cd_dens.rhotor_norm = abc.grids_1d.rhotornorm;
ohm_data.cd_dens.rhopol_norm = abc.grids_1d.rhopolnorm;
ohm_help = 'from IBS ohm related nodes, from Iohm=Ip-Icd-Ibs assuming stationary state';
else
disp(['source_icd.ohm = ' source_icd.ohm ' not yet implemented, ask O. Sauter'])
for subfields={'cd_tot','cd_dens'};
for subsubfields={'data','x','t','units','dim','dimunits','label','rhotor_norm','rhopol_norm'}
ohm_data.(subfields{1}).(subsubfields{1}) = [];
end
end
end
gdat_data.ohm.ohm_data = ohm_data;
data_fullpath = ['from ibs:ohmic related nodes, all results in .ohm_data, subfield=' field_for_main_data ...
'in ohm.data, .x, .t, .dim, .dimunits, .label, .units'];
for subfields={'data','x','t','units','dim','dimunits','label'}
gdat_data.ohm.(subfields{1}) = gdat_data.ohm.ohm_data.(field_for_main_data).(subfields{1});
end
if isempty(gdat_data.t); gdat_data.t = gdat_data.ohm.t; end
gdat_data.ohm.data_fullpath = data_fullpath;
gdat_data.ohm.help = ohm_help;
ij_ok = [isfinite(gdat_data.ohm.data)];
gdat_data.data(end+1,:) = interpos(-63,gdat_data.ohm.t(ij_ok),gdat_data.ohm.data(ij_ok),gdat_data.t,-0.1);
gdat_data.label{end+1}=gdat_data.ohm.label;
end
%
if any(strmatch('bs',gdat_data.gdat_params.source))
data_fullpath = '';
bs_help = '';
if strcmp(lower(source_icd.bs),'ibs')
bs_data.cd_tot = gdat_tcv(shot,'\results::ibs:ibs');
if gdat_data.mapping_for.tcv.gdat_timedim ==2
tgrid_to_change = {'bs_data.cd_tot.t','bs_data.cd_tot.dim{1}'};
for i=1:length(tgrid_to_change)
eval([tgrid_to_change{i} ' = reshape(' tgrid_to_change{i} ',1,numel(' tgrid_to_change{i} '));']);
end
end
bs_data.cd_tot.data = bs_data.cd_tot.data';
bs_data.cd_tot.units = strrep(bs_data.cd_tot.units,'kA','A');
bs_data.cd_dens = gdat_tcv(shot,'\results::ibs:jbsav'); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R
if gdat_data.mapping_for.tcv.gdat_timedim ==2
tgrid_to_change = {'bs_data.cd_dens.t','bs_data.cd_dens.dim{2}'};
for i=1:length(tgrid_to_change)
eval([tgrid_to_change{i} ' = reshape(' tgrid_to_change{i} ',1,numel(' tgrid_to_change{i} '));']);
end
end
bs_data.cd_dens.units = strrep(bs_data.cd_dens.units,'kA','A');
abc=get_grids_1d(bs_data.cd_dens,1,1);
bs_data.cd_dens.rhotor_norm = abc.grids_1d.rhotornorm;
bs_data.cd_dens.rhopol_norm = abc.grids_1d.rhopolnorm;
bs_help = 'from IBS bs related nodes, from Ibs=Ip-Icd-Ibs assuming stationary state';
else
disp(['source_icd.bs = ' source_icd.bs ' not yet implemented, ask O. Sauter'])
for subfields={'cd_tot','cd_dens'};
for subsubfields={'data','x','t','units','dim','dimunits','label','rhotor_norm','rhopol_norm'}
bs_data.(subfields{1}).(subsubfields{1}) = [];
end
end
end
gdat_data.bs.bs_data = bs_data;
data_fullpath = ['from ibs:bsic related nodes, all results in .bs_data, subfield=' field_for_main_data ...
'in bs.data, .x, .t, .dim, .dimunits, .label, .units'];
for subfields={'data','x','t','units','dim','dimunits','label'}
gdat_data.bs.(subfields{1}) = gdat_data.bs.bs_data.(field_for_main_data).(subfields{1});
end
if isempty(gdat_data.t); gdat_data.t = gdat_data.bs.t; end
gdat_data.bs.data_fullpath = data_fullpath;
gdat_data.bs.help = bs_help;
ij_ok = [isfinite(gdat_data.bs.data)];
gdat_data.data(end+1,:) = interpos(-63,gdat_data.bs.t(ij_ok),gdat_data.bs.data(ij_ok),gdat_data.t,-0.1);
gdat_data.label{end+1}=gdat_data.bs.label;
end
%
% add all to last index of .data(:,i)
gdat_data.data(end+1,:) = sum(gdat_data.data,1,'omitnan');
gdat_data.x = [1:size(gdat_data.data,1)];
gdat_data.label{end+1}='total';
gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits = {['index for each main source + total ' field_for_main_data], 's'};
gdat_data.data_fullpath = 'see in individual source substructure';
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
for i=1:length(gdat_data.gdat_params.source)
if ~isempty(gdat_data.(gdat_data.gdat_params.source{i}).data) && ~isempty(gdat_data.(gdat_data.gdat_params.source{i}).t)
% data with same time length as gdat_data.(gdat_data.gdat_params.source{i}).t
ab_tmp_rho = gdat_data.(gdat_data.gdat_params.source{i});
ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for;
switch gdat_data.gdat_params.source{i}
case 'ec'
other_data = {'ec_data.p_abs_plasma.data','ec_data.p_abs_plasma.t', ...
['ec_data.p_abs_plasma.dim{' num2str(gdat_data.mapping_for.tcv.gdat_timedim) '}'], ...
'ec_data.cd_tot.data','ec_data.cd_tot.t', ...
['ec_data.cd_tot.dim{' num2str(gdat_data.mapping_for.tcv.gdat_timedim) '}']};
case 'ohm'
other_data = {'ohm_data.cd_tot.data','ohm_data.cd_tot.t', ['ohm_data.cd_tot.dim{1}']};
case 'bs'
other_data = {'bs_data.cd_tot.data','bs_data.cd_tot.t', ['bs_data.cd_tot.dim{1}']};
otherwise
disp(gdat_data.gdat_params.source{i})
keyboard
end
[ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,other_data);
gdat_data.(gdat_data.gdat_params.source{i}) = rmfield(ab_tmp_rho,{'gdat_params','mapping_for'});
end
end
% extra ec subfields
if ~isempty(gdat_data.ec.data) && ~isempty(gdat_data.ec.t) && isfield(gdat_data.ec,'ec_data')
extra_fields=fieldnames(gdat_data.ec.ec_data);
extra_fields = setdiff(extra_fields,{'p_abs_plasma','cd_tot'});
for i=1:length(extra_fields)
ab_tmp_rho = gdat_data.ec.ec_data.(extra_fields{i});
ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for;
ab_tmp_rho.mapping_for.tcv.gdat_timedim = length(size(ab_tmp_rho.data)); % time always last dim for these fields
[ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1);
gdat_data.ec.ec_data.(extra_fields{i}) = rmfield(ab_tmp_rho,{'gdat_params','mapping_for'});
end
end
% extra cd_dens fields in ohm, bs
extra_fields = intersect({'bs','ohm'},gdat_data.gdat_params.source);
for i=1:length(extra_fields)
if ~isempty(gdat_data.(extra_fields{i}).data) && ~isempty(gdat_data.(extra_fields{i}).t) && ...
isfield(gdat_data.(extra_fields{i}),[extra_fields{i} '_data'])
gdat_data.(extra_fields{i}).([extra_fields{i} '_data']).cd_dens.gdat_params.time_out = gdat_data.gdat_params.time_out;
gdat_data.(extra_fields{i}).([extra_fields{i} '_data']).cd_dens = gdat2time_out( ...
gdat_data.(extra_fields{i}).([extra_fields{i} '_data']).cd_dens,1,{'rhotor_norm'});
end
end
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'mhd'}
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
warning(['parameter time_out not implemented yet, will load full time data; ask O. Sauter if needed' char(10)]);
gdat_data.gdat_params = rmfield(gdat_data.gdat_params,'time_out');
end
if isfield(gdat_data.gdat_params,'nfft') && ~isempty(gdat_data.gdat_params.nfft)
% used in gdat_plot for spectrogram plot
else
gdat_data.gdat_params.nfft = 1024;
end
% load n=1, 2 and 3 Bdot from magnetic measurements
n3.data = [];
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';
if shot>= 50926
gdat_data.dimunits{2} = 'n number, at this stage n3=n1';
gdat_data.dimunits{2} = 'n number, at this stage n3 not computed';
end
gdat_data.units = 'T/s';
gdat_data.request_description = 'delta_Bdot from magnetic probes to get n=1/odd, 2/even and 3';
gdat_data.label = {'n=1','n=2','n=3'}; % can be used in legend(gdat_data.label)
if shot>= 50926
gdat_data.label = {'n odd','n even'}; % can be used in legend(gdat_data.label)
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ne','te'}
% ne or Te from Thomson data on raw z mesh vs (z,t)
if ~(isfield(gdat_data.gdat_params,'systems') && ~isempty(gdat_data.gdat_params.systems) && ...
all(ismember(gdat_data.gdat_params.systems,{'all','main','edge'})))
gdat_data.gdat_params.systems = {'all'};
end
[gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.systems,gdat_params.nverbose,0);
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1,{'error_bar'});
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'nel'}
% line-averaged density
sources_avail = {'fir','fit'};
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = 'fir';
end
gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source);
if ~any(strmatch(gdat_data.gdat_params.source,sources_avail,'exact'))
error([gdat_data.gdat_params.source ' not defined, sources_avail = ' sources_avail]);
end
switch gdat_data.gdat_params.source
case 'fit'
psitbx01=psitbxtcv(gdat_data.shot,'01',psitbx_str);
neft = gdat_tcv(gdat_data.shot,'ne_rho','fit',1);
gdat_data.t = neft.fit.t;
gdat_data.dim{1} = gdat_data.t;
gdat_data.dimunits{1} = 's';
gdat_data.units = 'particles/m^3';
nel.dim{1} = neft.fit.t';
nel.data = ones(size(nel.dim{1}));
gdat_data.r_center_chord = 0.903;
[ratio,fir_times,z_lineint]=fir_ratio_from_fits(neft.fit.data,neft.fit.x,psitbx01,nel,gdat_data.r_center_chord,gdat_data.t',1,1e3);
gdat_data.data = reshape(1./ratio./z_lineint,numel(ratio),1);
gdat_data.data_fullpath = ['line-averaged density from fit using ' neft.fit.data_fullpath];
otherwise
gdat_data.gdat_params.data_request = '\results::fir:n_average';
gdat_data = gdat_tcv(gdat_data.shot,gdat_data.gdat_params);
ab_index = gdat_tcv(gdat_data.shot,'\diagz::fir_array:center_chord');
ab_r = gdat_tcv(gdat_data.shot,'\diagz::fir_array:radii');
gdat_data.r_center_chord = ab_r.data(ab_index.data+1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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,'systems') && ~isempty(gdat_data.gdat_params.systems) && ...
all(ismember(gdat_data.gdat_params.systems,{'all','main','edge'})))
gdat_data.gdat_params.systems = {'all'};
end
default_fit_type = 'conf';
if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit)
gdat_data.gdat_params.fit = 1;
end
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
data_option_for_firrat = 0; % note assume conf with local; should get from conf nodes the effective one
if strcmp(gdat_data.gdat_params.fit_type,'avg'); data_option_for_firrat = 1; end
[gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.systems,gdat_params.nverbose);
% construct rho mesh
% psiscatvol obtained from linear interpolation in fact so not quite ok near axis where psi is quadratic
% NOTE: Since 2019 psi_scatvol uses psitbx with cubic-spline interpolation, unsure what shots have new data or not ...
recompute_psiscatvol_always = 1;
if liuqe_version_eff==-1; recompute_psiscatvol_always = 1; end
% NOTE: Since oct.2019 psi_scatvol stored in nodes now uses LIUQE.M (shot range to check if any ...)
if all(abs(zshift)<1e-5) && liuqe_matlab==0 && recompute_psiscatvol_always==0
if strcmp(substr_liuqe,'_1'); substr_liuqe = ''; end
[psiscatvol,psi_max] = get_thomson_psiscatvol(shot,gdat_data.gdat_params.systems,substr_liuqe,gdat_params.nverbose);
% NOTE: time and Z are switched here (same as previous state) and
% code is not reachable anyway) ... waiting for a fix
else
% calculate new psiscatvol
[psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, psitbx_str, gdat_params);
end
% Add the results to the output of gdat
gdat_data.psiscatvol = psiscatvol;
gdat_data.psi_max = psi_max;
if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data)
psi_norm = 1.-psiscatvol.data./repmat(psi_max.data(:).',[size(psiscatvol.data,1),1]);
mask = psi_norm < 0;
if any(mask(:))
% NOTE: I am not entirely sure there could not be a case where this
% is actually possible outside of the confined plasma
warning('Found negative values of psi_norm at TS positions, replacing them with 0');
psi_norm(mask) = 0;
end
params_eff = gdat_data.gdat_params;
params_eff.doplot=0;
params_eff.data_request='r_contour_edge';
r_edge = gdat(gdat_data.shot,params_eff);
params_eff.data_request='z_contour_edge';
z_edge = gdat(gdat_data.shot,params_eff);
gdat_data.z_thomson = gdat_data.dim{1};
psinorm_max = max(max(psi_norm));
for it=1:numel(psiscatvol.dim{2})
iteq=iround_os(r_edge.t,psiscatvol.dim{2}(it));
inpol=inpolygon(gdat_data.r_thomson,gdat_data.z_thomson,r_edge.data(:,iteq),z_edge.data(:,iteq));
is_private_flux = [psi_norm(:,it) < 1] & ~inpol;
% at this stage set private flux "rho" to max, could use 1.3^2
psi_norm(is_private_flux,it) = psinorm_max;
end
rho = sqrt(psi_norm);
else
if gdat_params.nverbose>=1 && gdat_data.gdat_params.edge==0
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;
% add grids_1d to have rhotor, etc
gdat_data = get_grids_1d(gdat_data,2,1,gdat_params.nverbose);
%%%%%%%%%%%
% 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'});
%
te = get_thomson_raw_data(shot,'te',gdat_data,gdat_data.gdat_params.systems,gdat_params.nverbose);
gdat_data.te.data=te.data;
gdat_data.te.error_bar=te.error_bar;
gdat_data.te.units = te.units;
gdat_data.data_fullpath=['pe=1.6e-19*ne*Te in data, .ne, .te from \results::thomson for systems (' ...
strjoin(gdat_data.gdat_params.systems,'/') ') 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
% 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
switch gdat_data.gdat_params.fit_type
case 'avg'
def_proffit = '\results::proffit.avg_time:';
def_rhotornorm = '\results::proffit:rhotor';
def_rhovolnorm = '\results::proffit:rhovol';
case 'local'
def_proffit = '\results::proffit.local_time:';
def_rhotornorm = '\results::proffit:rhotor';
def_rhovolnorm = '\results::proffit:rhovol';
case 'conf'
def_proffit = '\results::conf:';
def_rhotornorm = '\results::conf:rhotor';
def_rhovolnorm = '\results::conf:rhovol';
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;
gdat_data.fit.help = '.x is rhopol coordinate';
tmp = tdi(def_rhotornorm);
gdat_data.fit.rhotornorm = tmp.data;
tmp = tdi(def_rhovolnorm);
gdat_data.fit.rhovolnorm = tmp.data;
% do te as well if nete asked for
if strcmp(data_request_eff(1:4),'nete')
for subfields={'data','x','t','units','help','rhotornorm','rhovolnorm','data_fullpath'}
gdat_data.fit.ne.(subfields{1}) = gdat_data.fit.(subfields{1});
end
% $$$ 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;
% $$$ gdat_data.fit.ne.help = gdat_data.fit.help;
% $$$ gdat_data.fit.ne.rhotornorm = gdat_data.fit.rhotornorm;
% $$$ gdat_data.fit.ne.rhovolnorm = gdat_data.fit.rhovolnorm;
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
for subfields={'x','t','help','rhotornorm','rhovolnorm'}
gdat_data.fit.te.(subfields{1}) = gdat_data.fit.ne.(subfields{1});
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
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
% add extra fields then correct
ab_tmp_rho = gdat_data.fit; ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for;
[ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,{'rhotornorm','rhovolnorm'});
gdat_data.fit = rmfield(ab_tmp_rho,{'gdat_params','mapping_for'});
%
gdat_data.psi_max.dim{2} = gdat_data.psi_max.dim{1};
extra_fields = {'error_bar','x','dim{1}','psiscatvol.data','psiscatvol.dim','psi_max.data','psi_max.dim'};
if strcmp(data_request_eff,'ne_rho')
extra_fields(end+1:end+3) = {'firrat','data_raw','error_bar_raw'};
gdat_data.firrat = reshape(gdat_data.firrat,1,numel(gdat_data.firrat));
end
[gdat_data] = gdat2time_out(gdat_data,22,extra_fields);
gdat_data.psi_max.dim = gdat_data.psi_max.dim{2};
if strcmp(data_request_eff,'nete_rho')
extra_sub = {'ne','te'};
for i=1:length(extra_sub)
ab_tmp_rho = gdat_data.fit.(extra_sub{i}); ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for;
[ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,{'rhotornorm','rhovolnorm'});
gdat_data.fit.(extra_sub{i}) = rmfield(ab_tmp_rho,{'gdat_params','mapping_for'});
%
ab_tmp_rho = gdat_data.(extra_sub{i});
ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for;
ab_tmp_rho.t = gdat_data.t;
extra_fields = {'error_bar'};
if strcmp(extra_sub{i},'ne');
extra_fields(end+1:end+3) = {'firrat','data_raw','error_bar_raw'};
ab_tmp_rho.firrat = reshape(ab_tmp_rho.firrat,1,numel(ab_tmp_rho.firrat));
end
[ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,extra_fields);
gdat_data.(extra_sub{i}) = rmfield(ab_tmp_rho,{'gdat_params','mapping_for','t'});
end
end
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'powers'}
% note: same time array for all main, ec, ohm, nbi, ...
% At this stage fill just ech, later add nbi
sources_avail = {'ohm','ec','nbi','rad'}; % note should allow ech, nbh, ohmic in parameter sources
% create empty structures for all, so in return one always have same substructres
for i=1:length(sources_avail)
gdat_data.(sources_avail{i}).data = [];
gdat_data.(sources_avail{i}).units = [];
gdat_data.(sources_avail{i}).dim=[];
gdat_data.(sources_avail{i}).dimunits=[];
gdat_data.(sources_avail{i}).t=[];
gdat_data.(sources_avail{i}).x=[];
gdat_data.(sources_avail{i}).data_fullpath=[];
gdat_data.(sources_avail{i}).label=[];
end
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = sources_avail;
elseif ~iscell(gdat_data.gdat_params.source)
if ischar(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source);
if ~any(strmatch(gdat_data.gdat_params.source,lower(sources_avail)))
if (gdat_params.nverbose>=1)
warning(['source= ' gdat_data.gdat_params.source ' is not part of the available sources: ' sprintf('''%s'' ',sources_avail{:})]);
end
return
else
gdat_data.gdat_params.source = {gdat_data.gdat_params.source};
end
else
if (gdat_params.nverbose>=1); warning([' source parameter not compatible with: ' sprintf('''%s'' ',sources_avail{:})]); end
return
end
else
for i=1:length(gdat_data.gdat_params.source)
gdat_data.gdat_params.source{i} = lower(gdat_data.gdat_params.source{i});
if ~any(strmatch(gdat_data.gdat_params.source{i},lower(sources_avail)))
if gdat_data.gdat_params.nverbose>=1
warning(['source = ' gdat_data.gdat_params.source{i} ' not expected with data_request= ' data_request_eff])
end
end
end
end
% always start from ohmic so can use this time as base time since should yield full shot
% get ohmic power simply from vloop*Ip (minus sign for TCV), neglect dWp/dt could add it later, see chie_tcv_to_nodes
% thus should take it from conf if present
mdsopen(shot);
ptot_ohm = tdi('\results::conf:ptot_ohm');
if ~isempty(ptot_ohm.data) && ~ischar(ptot_ohm.data) && ~isempty(ptot_ohm.dim)
gdat_data.ohm.data = ptot_ohm.data;
gdat_data.ohm.t = ptot_ohm.dim{1};
gdat_data.ohm.dim = ptot_ohm.dim;
gdat_data.ohm.dimunits = ptot_ohm.dimunits;
gdat_data.ohm.units = ptot_ohm.units;
gdat_data.ohm.data_fullpath = '\results::conf:ptot_ohm';
gdat_data.ohm.help = ptot_ohm.help;
else
params_eff = gdat_data.gdat_params;
params_eff.data_request='ip'; % to make sure to use input params like liuqe option
ip=gdat_tcv([],params_eff); %gdat_tcv to avoid plotting in case doplot=1 if using gdat and to save time
params_eff.data_request='vloop';
vloop=gdat_tcv([],params_eff);
gdat_data.ohm.t = vloop.t;
gdat_data.ohm.dim{1} = gdat_data.t;
gdat_data.dimunits{1} = 's';
gdat_data.units = 'W';
tension = -1e5;
vloop_smooth=interpos(-63,vloop.t,vloop.data,gdat_data.ohm.t,tension);
ip_t = interp1(ip.t,ip.data,gdat_data.ohm.t);
gdat_data.ohm.data = -vloop_smooth.*ip_t; % TCV has wrong sign for Vloop
gdat_data.ohm.raw_data = -vloop.data.*ip_t;
gdat_data.ohm.data_fullpath = 'from vloop*Ip, smoothed vloop in data, unsmoothed in raw_data';
gdat_data.ohm.data_fullpath=['from -vloop(tens=' num2str(tension,'%.0e') ')*ip, from gdat, in .data, unsmoothed in .raw_data'];
end
gdat_data.ohm.x=[];
gdat_data.ohm.label='P_{ohm}';
%
% add each source in main.data, on ohm time array
gdat_data.t = linspace(gdat_data.ohm.t(1),gdat_data.ohm.t(end),floor(1e4.*(gdat_data.ohm.t(end)-gdat_data.ohm.t(1))))';
ij=[isfinite(gdat_data.ohm.data)];
gdat_data.data(:,1) = interpos(-21,gdat_data.ohm.t(ij),gdat_data.ohm.data(ij),gdat_data.t);
gdat_data.dim{1} = gdat_data.t;
gdat_data.x(1) = 1;
gdat_data.label={'P_{ohm}'};
gdat_data.units = 'W';
%
if any(strmatch('ec',gdat_data.gdat_params.source))
% EC
nodenameeff='\results::toray.input:p_gyro';
tracetdi=tdi(nodenameeff);
if isempty(tracetdi.data) || isempty(tracetdi.dim) || ischar(tracetdi.data)
if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
else
ij=[~isfinite(tracetdi.data)]; tracetdi.data(ij)=0.;
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}';
gdat_data.ec.help = tracetdi.help;
% add to main with linear interpolation and 0 for extrapolated values
gdat_data.data(:,end+1) = interpos(-21,gdat_data.ec.t,gdat_data.ec.data(:,end),gdat_data.t);
gdat_data.x(end+1) = size(gdat_data.data,2);
gdat_data.label{end+1}=gdat_data.ec.label;
end
end
%
if any(strmatch('nb',gdat_data.gdat_params.source))
NBH_in_TCV = 0;
if shot >= 51641
% NBI
moderemote = mdsdata('data(\VSYSTEM::TCV_NBHDB_B["NBI:MODEREMOTE_FB"])');
lcs_mode = mdsdata('data(\VSYSTEM::TCV_NBHDB_I["NBI:LCS_MODE"])');
% NBH in TCV equiv moderemote='ON' AND lcs_mode = 9
NBH_in_TCV = strcmpi(strtrim(moderemote),'ON') && lcs_mode == 9;
else
% Nodes used in previous block only exist outside of Vista for shots after 51641
if any(shot == [51458 51459 51460 51461 51463 51465 51470 51472 ... % 29.JAN.2016
51628 51629 51631 51632 51633 ... % 09.FEB.2016
51639 51640 ... 51641 % 10.FEB.2016
])
NBH_in_TCV = 1;
end
end
if NBH_in_TCV
nodenameeff = '\results::NBH:POWR_TCV';
nbh_data_tdi = tdi(nodenameeff);
if ~isempty(nbh_data_tdi.data) && ~ischar(nbh_data_tdi.data) && ~isempty(nbh_data_tdi.dim)
nbi_neutral_power_tot = nbh_data_tdi.data.*1e6; % in W
ij = nbi_neutral_power_tot<100;
nbi_neutral_power_tot(ij) = 0.;
gdat_data.nbi.data = nbi_neutral_power_tot; % at this stage p_gyro is in kW'
gdat_data.nbi.units = 'W';
gdat_data.nbi.dim{1}=nbh_data_tdi.dim{1};
gdat_data.nbi.dimunits{1}=nbh_data_tdi.dimunits{1};
gdat_data.nbi.t=gdat_data.nbi.dim{1};
gdat_data.nbi.x=[];
gdat_data.nbi.data_fullpath=[nodenameeff];
gdat_data.nbi.label='P_{nbi}';
gdat_data.nbi.help = nbh_data_tdi.help;
nodenameeff2 = '\results::nbh:energy';
gdat_data.nbi.data_fullpath=[nodenameeff ' *1e6 for W; and ' nodenameeff2 ' *1e3 for eV'];
nbh_energy_data_tdi = tdi(nodenameeff2);
gdat_data.nbi.energy = nbh_energy_data_tdi.data*1e3; % for eV
% add to main with linear interpolation and 0 for extrapolated values
ij=[isfinite(gdat_data.nbi.data)];
gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi.t(ij),gdat_data.nbi.data(ij),gdat_data.t);
gdat_data.x(end+1) = size(gdat_data.data,2);
gdat_data.label{end+1}=gdat_data.nbi.label;
end
end
end
%
index_rad = [];
if any(strmatch('rad',gdat_data.gdat_params.source))
% RAD
nodenameeff='\results::bolo:prad:total';
tracetdi=tdi(nodenameeff);
if isempty(tracetdi.data) || isempty(tracetdi.dim) || ischar(tracetdi.data)
if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
else
gdat_data.rad.data = tracetdi.data*1e3; % at this stage bolo is in kW'
gdat_data.rad.units = 'W';
gdat_data.rad.dim=tracetdi.dim;
gdat_data.rad.dimunits=tracetdi.dimunits;
gdat_data.rad.t=tracetdi.dim{1};
gdat_data.rad.data_fullpath=[nodenameeff];
gdat_data.rad.label='P_{RAD}';
gdat_data.rad.help = tracetdi.help;
% add to main with linear interpolation and 0 for extrapolated values
gdat_data.data(:,end+1) = interpos(-21,gdat_data.rad.t,gdat_data.rad.data(:,end),gdat_data.t);
index_rad = size(gdat_data.data,2);
gdat_data.x(end+1) = size(gdat_data.data,2);
gdat_data.label{end+1}=gdat_data.rad.label;
end
end
% add all to last index of .data(:,i)
ij = setdiff([1:size(gdat_data.data,2)],index_rad);
gdat_data.data(:,end+1) = sum(gdat_data.data(:,ij),2);
gdat_data.x(end+1) = size(gdat_data.data,2);
gdat_data.label{end+1}='total heating';
gdat_data.dim{2} = gdat_data.x;
gdat_data.dimunits = {'s', 'index for each source + total heating'};
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
for i=1:length(sources_avail)
ab_tmp_rho = gdat_data.(sources_avail{i});
ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for;
if strcmp(sources_avail{i},'nbi')
[ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,{'energy'});
else
[ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1);
end
gdat_data.(sources_avail{i}) = rmfield(ab_tmp_rho,{'gdat_params','mapping_for'});
end
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'q_rho'}
% q profile on psi from liuqe
if liuqe_matlab==0
nodenameeff=['tcv_eq(''q_psi'',''LIUQE' substr_liuqe_tcv_eq ''')'];
else
nodenameeff=['tcv_eq(''q'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
end
if liuqe_version_eff==-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'];
if liuqe_matlab==0
rhopol_eff = ones(size(tracetdi.dim{1}));
rhopol_eff(:) = sqrt(linspace(0,1,length(tracetdi.dim{1})));
gdat_data.dim{1} = rhopol_eff;
end
gdat_data.x = gdat_data.dim{1};
gdat_data.dimunits{1} = 'rho_pol~sqrt(\psi_norm)';
gdat_data.dimunits{2} = 's';
gdat_data.units = '';
gdat_data.request_description = 'q(rhopol\_norm)';
% add grids_1d to have rhotor, etc
gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose);
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,21);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'rbphi_rho', 'rbtor_rho'}
% R*Bphi(rho,t) from F from FFprime
if liuqe_matlab==0
disp('not yet implemented for liuqe fortran')
return
else
nodenameeff=['tcv_eq(''rbtor_rho'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
end
if liuqe_version_eff==-1
disp('not yet implemented for fbte')
return
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'];
if liuqe_matlab==0
rhopol_eff = ones(size(tracetdi.dim{1}));
rhopol_eff(:) = sqrt(linspace(0,1,length(tracetdi.dim{1})));
gdat_data.dim{1} = rhopol_eff;
end
gdat_data.x = gdat_data.dim{1};
gdat_data.dimunits{1} = 'rho_pol~sqrt(\psi_norm)';
gdat_data.dimunits{2} = 's';
gdat_data.units = '';
gdat_data.request_description = nodenameeff;
% add grids_1d to have rhotor, etc
gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose);
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
% add extra fields then correct
[gdat_data] = gdat2time_out(gdat_data,21);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
case {'phi_tor', 'phitor', 'toroidal_flux'}
% 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
% need to get q_rho but to avoid loop for rhotor in grids_1d, get q_rho explicitely here
params_eff = gdat_data.gdat_params;
if liuqe_matlab==1
nodenameeff = ['tcv_eq(''tor_flux_tot'',''' psitbx_str ''')'];
phi_tor = tdi(nodenameeff);
else
phi_tor.data = [];
phi_tor.units = 'Wb';
end
if ~any(any(isfinite(phi_tor.data))) || ischar(phi_tor.data)
% no phi_tor, compute it from q profile
q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']);
if liuqe_matlab==0
q_rho.dim{1} = sqrt(linspace(0,1,numel(q_rho.dim{1})));
end
phi_tor.dim = q_rho.dim;
phi_tor.dimunits = q_rho.dimunits;
params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE
psi_axis=gdat_tcv(shot,params_eff);
if isnumeric(q_rho.data)
for it=1:size(q_rho.data,2)
ij=find(isfinite(q_rho.data(:,it)));
if ~isempty(ij) && numel(ij)>5
[qfit,~,~,phi]=interpos(-43,q_rho.dim{1}(ij).^2,q_rho.data(ij,it),phi_tor.dim{1}.^2);
% Phi=sigma_bp*sigma_rhothetaphi*(2pi)^(1-eBp)*int(q dpsi), dpsi=dpsi_norm * (psiedge-psaxis)
% cocos=17: sigma_Bp=-1,sigma_rhothetaphi=1, eBp=1, (psiedge-psaxis)=-psiaxis
phi = -1.* phi .* (-psi_axis.data(it));
phi_tor.data(:,it) = phi';
end
end
else
phi_tor.data = [];
end
elseif any(any(~isfinite(phi_tor.data)))
% there are some non-finite values, usually edge values, so replace them with integrated values
only_edge = 0;
if any(any(~isfinite(phi_tor.data(1:end-1,:))))
only_edge = 1;
end
q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']);
if liuqe_matlab==0
q_rho.dim{1} = sqrt(linspace(0,1,numel(q_rho.dim{1})));
end
params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE
psi_axis=gdat_tcv(shot,params_eff);
% assume same dimensions, check. Note data dims can be different from dim{} when there is a problem with liuqe run
if numel(phi_tor.data) ~= numel(q_rho.data) || numel(psi_axis.data) ~= numel(phi_tor.dim{2})
warning(['problems in gdat_tcv with ' data_request_eff ' with unexpected dimensions'])
return;
end
% check time sizes since can happen problems with old liuqe runs
if size(q_rho.data,2) ~= size(phi_tor.data,2)
warning(['time sizes between q_rho and phi_tor do not match: ' num2str(size(q_rho.data,2)) ' and ' num2str(size(phi_tor.data,2))]);
return
end
if size(q_rho.data,2) ~= length(psi_axis.data)
psi_axis
q_rho
disp('q_rho.dim');q_rho.dim
phi_tor
disp('phi_tor.dim');phi_tor.dim
disp(['WARNING: time sizes between q_rho and psi_axis do not match: ' num2str(size(q_rho.data,2)) ' and ' ...
num2str(length(psi_axis.data))]);
disp('WARNING: should have been caught, so probably q_rho.dim not the same as data(:,:) dims')
disp('should run: >> liuqe_default_run_sequence(shot)')
return
end
if size(q_rho.data,1) ~= size(phi_tor.data,1)
warning(['radial sizes between q_rho and phi_tor do not match: ' num2str(size(q_rho.data,1)) ' and ' num2str(size(phi_tor.data,1))]);
return
end
for it=1:size(q_rho.data,2)
%if any(~isfinite(phi_tor.data(:,it)))
ij=find(isfinite(q_rho.data(:,it)));
if ~isempty(ij) && numel(ij)>5
[qfit,~,~,phi]=interpos(-43,q_rho.dim{1}(ij).^2,q_rho.data(ij,it),phi_tor.dim{1}.^2);
phi = -1.* phi .* (-psi_axis.data(it));
if only_edge
phi_tor.data(end,it) = phi(end);
else
phi_tor.data(:,it) = phi';
end
end
%end
end
end
gdat_data.data = phi_tor.data;
gdat_data.dim = phi_tor.dim;
gdat_data.dimunits = phi_tor.dimunits;
gdat_data.units = phi_tor.units;
if (length(gdat_data.dim)>=mapping_for_tcv.gdat_timedim); gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim}; end
if (length(gdat_data.dim)>0); gdat_data.x = gdat_data.dim{1}; end
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
% add extra fields then correct
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'pprime', 'pprime_rho'}
if liuqe_matlab==0
nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('pprime_rho',1,0) '",''' psitbx_str ''')'];
tracetdi=tdi(nodenameeff);
if ~isempty(tracetdi.dim) && length(tracetdi.dim)>=1
tracetdi.dim{1} = sqrt(tracetdi.dim{1}); % correct x-axis psi_norm to rhopol
end
tracetdi.data = tracetdi.data ./2 ./pi; % correct node assumption (same for liuqe_fortran and fbte)
else
nodenameeff = ['tcv_eq(''pprime_rho'',''' psitbx_str ''')'];
tracetdi=tdi(nodenameeff);
end
gdat_data.data = tracetdi.data;
gdat_data.dim = tracetdi.dim;
if ~isempty(gdat_data.dim) && length(gdat_data.dim)>=2
gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
gdat_data.x = gdat_data.dim{setdiff([1 2],mapping_for_tcv.gdat_timedim)};
end
gdat_data.data_fullpath=nodenameeff;
gdat_data.dimunits = tracetdi.dimunits;
gdat_data.units = tracetdi.units;
gdat_data.request_description = 'pprime=dp/dpsi';
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
case {'pressure', 'pressure_rho', 'p_rho'}
if liuqe_matlab==0
nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('p_rho',1,0) '",''' psitbx_str ''')'];
tracetdi=tdi(nodenameeff);
if ~isempty(tracetdi.dim) && length(tracetdi.dim)>=1
tracetdi.dim{1} = sqrt(tracetdi.dim{1}); % correct x-axis psi_norm to rhopol
end
tracetdi.data = tracetdi.data ./2 ./pi; % correct node assumption (same for liuqe_fortran and fbte)
else
nodenameeff = ['tcv_eq(''p_rho'',''' psitbx_str ''')'];
tracetdi=tdi(nodenameeff);
end
gdat_data.data = tracetdi.data;
gdat_data.dim = tracetdi.dim;
if ~isempty(gdat_data.dim) && length(gdat_data.dim)>=2
gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
gdat_data.x = gdat_data.dim{setdiff([1 2],mapping_for_tcv.gdat_timedim)};
end
gdat_data.data_fullpath=nodenameeff;
gdat_data.dimunits = tracetdi.dimunits;
gdat_data.units = tracetdi.units;
gdat_data.request_description = 'pressure';
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'psi_edge'}
% psi at edge, 0 by construction in Liuqe, thus not given
% add surface_psi from surface_flux and d(surface_flux)/dt = vloop
warning('now psi_edge expression in tcv_requests_mapping, so should not come here')
keyboard
nodenameeff=['\results::psi_axis' substr_liuqe];
if liuqe_version_eff==-1
nodenameeff=[begstr 'q_psi' substr_liuqe];
end
tracetdi=tdi(nodenameeff);
if isempty(tracetdi.data) || isempty(tracetdi.dim) || ~any(isfinite(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';
%
params_eff.data_request='\results::surface_flux';
surface_psi=gdat_tcv([],params_eff);
ij=[isfinite(surface_psi.data)];
[aa,vsurf] = interpos(63,surface_psi.t(ij),surface_psi.data(ij),-3);
gdat_data.surface_psi = surface_psi.data;
gdat_data.vsurf = vsurf;
gdat_data.vsurf_description = ['time derivative from surface_psi, obtained from \results::surface_flux'];
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'r_contour_edge', 'z_contour_edge'}
if liuqe_matlab==0
nodenameeff=['\results::' data_request_eff(1) '_contour' substr_liuqe];
% npts_contour = tdi(['\results::npts_contour' substr_liuqe]);
tracetdi=tdi(nodenameeff);
gdat_data.request_description = 'NaNs when smaller nb of boundary points at given time, can use \results::npts_contour';
else
nodenameeff = ['tcv_eq(''' data_request_eff(1) '_edge'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
tracetdi=tdi(nodenameeff);
gdat_data.request_description = 'lcfs on same nb of points for all times';
end
gdat_data.data = tracetdi.data;
gdat_data.dim = tracetdi.dim;
gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
gdat_data.data_fullpath=nodenameeff;
gdat_data.dimunits = tracetdi.dimunits;
gdat_data.units = tracetdi.units;
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'rhotor_edge', 'rhotor', 'rhotor_norm'}
% 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
% need to get q_rho but to avoid loop for rhotor in grids_1d, get q_rho explicitely here
params_eff = gdat_data.gdat_params;
if liuqe_matlab==0
nodenameeff=['tcv_eq(''q_psi'',''LIUQE' substr_liuqe_tcv_eq ''')'];
if liuqe_version_eff==-1
nodenameeff=[begstr 'q_psi' substr_liuqe];
end
q_rho=tdi(nodenameeff);
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
rhopol_eff = ones(size(q_rho.dim{1}));
rhopol_eff(:) = sqrt(linspace(0,1,length(q_rho.dim{1})));
q_rho.dim{1} = rhopol_eff;
params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE
psi_axis=gdat_tcv(shot,params_eff);
params_eff.data_request='b0';
b0=gdat_tcv(shot,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 isfield(params_eff,'rhopol') && ~isempty(params_eff.rhopol)
if numel(params_eff.rhopol) ~= 1
rhoequal = params_eff.rhopol;
elseif params_eff.rhopol > 0
rhoequal = linspace(0,1,params_eff.rhopol); % equidistant rhopol with given nb points
elseif params_eff.rhopol < 0
rhoequal = sqrt(linspace(0,1,abs(params_eff.rhopol))); % equidistant psinorm with given nb points
end
rhoequal = reshape(rhoequal,1,numel(rhoequal));
else
params_eff.rhopol = []; % means keep default rho mesh
end
if strcmp(data_request_eff,'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_eff,'rhotor_norm') || strcmp(data_request_eff,'rhotor')
% keep orientation as original
gdat_data.dim{1} = shiftdim(rhoequal,abs(find(size(q_rho.dim{1})==1)-find(size(rhoequal)==1)));
gdat_data.dim{2} = q_rho.dim{2};
gdat_data.data = NaN(numel(gdat_data.dim{1}),numel(gdat_data.dim{2})); % to have the dimensions correct
gdat_data.t = gdat_data.dim{2};
if strcmp(data_request_eff,'rhotor_norm')
gdat_data.data_fullpath='sqrt(phitor/phitor_edge), rhotor_edge=sqrt(phitor(1)/B0/pi)';
else
gdat_data.data_fullpath='sqrt(phitor/pi/B0), rhotor_edge=sqrt(phitor(1)/B0/pi)';
end
gdat_data.units = '';
gdat_data.dimunits{1} = 'rhopol\_norm';
gdat_data.dimunits{2} = 's';
end
gdat_data.x=gdat_data.dim{1};
gdat_data.psi_axis = reshape(psi_axis.data,length(psi_axis.data),1);
if length(gdat_data.psi_axis) ~= length(gdat_data.t)
disp('problems of time between qrho and psi_axis?')
return
end
gdat_data.b0 = b0tpsi;
for it=1:size(q_rho.data,2)
ij=find(isfinite(q_rho.data(:,it)));
if ~isempty(ij)
[qfit,~,~,phi]=interpos(-43,q_rho.dim{1}(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
gdat_data.rhotor_edge(it) = dataeff(end); % redundant with "rhotor_edge" but to always have all subfields
if strcmp(data_request_eff,'rhotor_edge')
gdat_data.data(it) = dataeff(end);
elseif strcmp(data_request_eff,'rhotor_norm')
gdat_data.data(:,it) = dataeff./dataeff(end);
elseif strcmp(data_request_eff,'rhotor')
gdat_data.data(:,it) = dataeff;
else
disp(['problem in gdat_tcv rhotor with unknown data_request_eff = ' data_request_eff]);
return
end
end
gdat_data.rhotor_edge = gdat_data.rhotor_edge';
else
params_eff = gdat_data.gdat_params;
params_eff.data_request='phi_tor';
phi_tor=gdat_tcv([],params_eff);
params_eff = gdat_data.gdat_params;
params_eff.data_request='b0';
b0=gdat_tcv([],params_eff);
gdat_data.t = phi_tor.t;
ij=find(isfinite(b0.data));
if length(gdat_data.t) > 0
gdat_data.b0 = interpos(b0.t(ij),b0.data(ij),gdat_data.t);
% always provide rhotor_edge so field always exists
gdat_data.rhotor_edge = sqrt(phi_tor.data(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)))';
if strcmp(data_request_eff,'rhotor_edge')
gdat_data.data = gdat_data.rhotor_edge;
gdat_data.dim = phi_tor.dim(2); % while there is a problem with \tcv_shot::top.results.equil... dim nodes
gdat_data.dimunits = phi_tor.dimunits(2);
gdat_data.units = 'm';
gdat_data.data_fullpath=['rhotor_edge=sqrt(Phi(edge)/pi/B0) from phi_tor'];
elseif strcmp(data_request_eff,'rhotor_norm')
gdat_data.data = sqrt(phi_tor.data./(ones(size(phi_tor.data,1),1)*reshape(phi_tor.data(end,:),1,size(phi_tor.data,2))));
gdat_data.dim = phi_tor.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes
gdat_data.x=gdat_data.dim{1};
gdat_data.units = '';
gdat_data.dimunits = phi_tor.dimunits;
gdat_data.data_fullpath=['rhotor_norm=sqrt(Phi/Phi(edge)) from phi_tor'];
elseif strcmp(data_request_eff,'rhotor')
gdat_data.data = sqrt(phi_tor.data./pi./(ones(size(phi_tor.data,1),1)*reshape(gdat_data.b0,1,numel(gdat_data.b0))));
gdat_data.dim = phi_tor.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes
gdat_data.x=gdat_data.dim{1};
gdat_data.units = 'm';
gdat_data.dimunits = phi_tor.dimunits;
gdat_data.data_fullpath=['rhotor=sqrt(Phi/pi/B0) from phi_tor'];
else
disp(['data_request_eff = ' data_request_eff ' not defined within rhotor block in gdat_tcv.m']);
return
end
else
gdat_data.b0 = [];
gdat_data.rhotor_edge = [];
end
end
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
% nothing to do since above uses gdat call getting reduced time already
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'rhovol','rho_vol','volume_rho','volume'}
if strcmp(data_request_eff,'rho_vol'); data_request_eff = 'rhovol'; end
% volume_rho = vol(rho); volume = vol(LCFS) = vol(rho=1);
% rhovol = sqrt(vol(rho)/vol(rho=1));
if liuqe_matlab==0
nodenameeff=['\results::psitbx:vol'];
if liuqe_version_eff > 1
disp('needs to construct volume');
return
end
else
nodenameeff=['tcv_eq(''vol'',''' psitbx_str ''')'];
end
if liuqe_version_eff==-1
data_request_eff = 'volume'; % only LCFS
nodenameeff=[begstr 'volume' substr_liuqe];
end
tracetdi=tdi(nodenameeff);
if (isempty(tracetdi.data) || isempty(tracetdi.dim)) && liuqe_matlab==0
% try to run psitbxput
psitbxput_version = 1.3;
psitbxput(psitbxput_version,shot);
ishot = mdsopen(shot);
tracetdi=tdi(nodenameeff);
end
if isempty(tracetdi.data) || isempty(tracetdi.dim) || ischar(tracetdi.data)
return
end
gdat_data.units = tracetdi.units;
if strcmp(data_request_eff,'volume')
if length(tracetdi.dim) >= 2
gdat_data.data = reshape(tracetdi.data(end,:),numel(tracetdi.data(end,:)),1); % needed consistent with timedim for gdat2time
gdat_data.volume_edge = gdat_data.data;
gdat_data.dim{1} = reshape(tracetdi.dim{2},numel(tracetdi.dim{2}),1);
gdat_data.dimunits{1} = tracetdi.dimunits{2};
else
mapping_for_tcv.gdat_timedim = 1;
gdat_data.data = reshape(tracetdi.data,numel(tracetdi.data),1);
gdat_data.dim = reshape(tracetdi.dim,numel(tracetdi.dim),1);
gdat_data.dimunits = tracetdi.dimunits;
gdat_data.volume_edge = gdat_data.data;
end
gdat_data.request_description = 'volume(LCFS)=volume(rhopol=1)';
gdat_data.data_fullpath=['Volume of LCFS from ' nodenameeff];
else
gdat_data.data = tracetdi.data;
gdat_data.dim = tracetdi.dim;
if length(gdat_data.dim)>0; gdat_data.x = gdat_data.dim{1}; end
gdat_data.dimunits = tracetdi.dimunits;
% to always have field volume_edge
gdat_data.volume_edge = gdat_data.data(end,:);
if strcmp(data_request_eff,'volume_rho')
gdat_data.request_description = 'volume(rho)';
gdat_data.data_fullpath=['Volume(rho) from ' nodenameeff];
elseif strcmp(data_request_eff,'rhovol')
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(volume/volume(edge) from ' nodenameeff];
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};
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1,{'volume_edge'});
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'rtc'}
% load all real-time memory signals for various nodes
%Get the data from mds and fill the data structure defined by
%define_simulink_signals
% If the rtccode/tools not in paths use sauter's link
% need: get_scd_mems.m
%
aaa=which('get_scd_mems');
if isempty(aaa)
if exist('/home/sauter/NoTivoli/rtccode_develop/tools')
addpath('/home/sauter/NoTivoli/rtccode_develop/tools','-end');
else
warning('Cannot add path for get_scd_mems')
return
end
end
aa = get_scd_mems(shot);
if isstruct(aa)
scd_names = fieldnames(aa);
for i=1:length(scd_names)
gdat_data.scd_mems.(scd_names{i}) = aa.(scd_names{i});
end
end
mdsclose;
mdsdisconnect;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
time_interval = [];
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
if length(gdat_data.gdat_params.time_out) == 2
time_interval = gdat_data.gdat_params.time_out;
else
if length(gdat_data.gdat_params.time_out) == 1
% 200ms includes all characteristic time constants
time_interval = [gdat_data.gdat_params.time_out-0.1 gdat_data.gdat_params.time_out+0.1];
else
time_interval = [min(gdat_data.gdat_params.time_out)-0.1 max(gdat_data.gdat_params.time_out)+0.1];
end
warning(['Expects a time interval [t1 t2] for ' data_request_eff ' in param time_out, uses [' ...
num2str(time_interval(1)) ',' num2str(time_interval(2)) ']' char(10)])
end
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_eff==[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(time_interval); t_int = 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_eff,'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(time_interval);
[sig,t,Sat_channel,cat_default] = XTOMOGetData(shot,0,4,camera_xtomo,[]);
else
[sig,t,Sat_channel,cat_default] = XTOMOGetData(shot,time_interval(1),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 {'transp'}
% read transp netcdf output file provided in source parameter and generate substructure with all variables
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || ~exist(gdat_data.gdat_params.source,'file')
warning([char(10) 'requires a netcdf file in ''source'' parameter' char(10)])
return
end
gdat_data.transp = cdf2mat(gdat_data.gdat_params.source);
% deal with errors and exceptions
% XB and X are rhotor_norm and not r/a radii
vars2change = {'X', 'XB'};
main_substruct = {'coords','allvars'};
subfields = {'long_name', 'label'};
for i=1:length(vars2change)
for j=1:length(main_substruct)
for k=1:length(subfields)
try
gdat_data.transp.(main_substruct{j}).(vars2change{i}).(subfields{k}) = strrep(gdat_data.transp.coords.(vars2change{i}).(subfields{k}),'x"r/a" bdy','sqrt(Phi/Philim)');
gdat_data.transp.(main_substruct{j}).(vars2change{i}).(subfields{k}) = strrep(gdat_data.transp.coords.(vars2change{i}).(subfields{k}),'x"r/a" ctr', ...
'sqrt(Phi/Philim)');
catch
end
end
end
end
gdat_data.data_fullpath = ['.transp structure obtained from cdf2mat(gdat_data.gdat_params.source)'];
gdat_data.help='Can use transp_gui(gdat_data.transp) to browse and plot';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ttprime', 'ttprime_rho'}
if liuqe_matlab==0
nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('ttprime_rho',1,0) '",''' psitbx_str ''')'];
tracetdi=tdi(nodenameeff);
if ~isempty(tracetdi.dim) && length(tracetdi.dim)>=1
tracetdi.dim{1} = sqrt(tracetdi.dim{1}); % correct x-axis psi_norm to rhopol
end
tracetdi.data = tracetdi.data ./2 ./pi; % correct node assumption (same for liuqe_fortran and fbte)
else
nodenameeff = ['tcv_eq(''ttprime_rho'',''' psitbx_str ''')'];
tracetdi=tdi(nodenameeff);
end
gdat_data.data = tracetdi.data;
gdat_data.dim = tracetdi.dim;
if ~isempty(gdat_data.dim) && length(gdat_data.dim)>=2
gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
gdat_data.x = gdat_data.dim{setdiff([1 2],mapping_for_tcv.gdat_timedim)};
end
gdat_data.data_fullpath=nodenameeff;
gdat_data.dimunits = tracetdi.dimunits;
gdat_data.units = tracetdi.units;
gdat_data.request_description = 'ttprime=t*dt/dpsi';
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
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];
gdat_data.mapping_for.tcv.gdat_timedim = 2;
if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
[gdat_data] = gdat2time_out(gdat_data,1);
elseif ~isfield(gdat_data.gdat_params,'time_out')
gdat_data.gdat_params.time_out = [];
end
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,timebase,nverbose,data_option)
%
% return fir_to_thomson ratio on time=timebase
%
% Use of autofit or thomson fir ratio depends on shot number
%
% Update on May 5th 2020 (amerle):
% Node for edge system thomson.edge:fir_thom_rat was never used so this
% option is now removed. The edge system existed for shots 24382:49895
% which lies entirely in the range for the autofit ratio. Note also that a
% separate fir ratio for edge and main makes very little sense as none of
% them cover the whole chord height.
%
% Update on Aug 8 2020
% Should take fir ratio depending on data taken:
% From auto fit data, use fir_thom_rat
% Otherwise in general use from local_time proffit: \results::proffit.local_time:neft_firat
% If using avg_time data, use \results::proffit.avg_time:neft_firat
%
% data_option: 0 (or empty), default: \results::proffit.local_time:neft_firat (best fit local data)
% 1: from avg_time fit: \results::proffit.avg_time:neft_firat (if this is best fit or CONF uses avg_time)
% 2: from auto_fits: fir_thom_rat (old, deprecated)
firthomratio = NaN;
if ~exist('timebase','var') || isempty(timebase)
if (nverbose>=1)
disp('need a timebase in get_fir_thom_rat_data')
end
return
end
firthomratio = NaN*ones(size(timebase));
data_option_eff = 0;
if exist('data_option') && ~isempty(data_option) && isfinite(data_option); data_option_eff = data_option; end
switch data_option_eff
case 0
tracefirrat=tdi('\results::proffit.local_time:neft_firat');
case 1
tracefirrat=tdi('\results::proffit.avg_time:neft_firat');
case 2
warning('Should use ratio computed with ''best'' fit: either proffit:local_time (default) or proffit_avg_time (if CONF used avg_time)')
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
otherwise
error(['case data_option_eff = ' num2str(data_option_eff) ' not defined in gdat_tcv: get_fir_thom_rat_data'])
end
if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) && any(isfinite(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,systems,nverbose,data_option_for_firrat)
%
% Get ne and Te from Thomson
%
% Also gets fir to Thomson ratio, which should be based on best fit mapping:
% data_option_for_firrat: see get_fir_thom_rat_data for details and default
%
if ~exist('data_option_for_firrat')
data_option_for_firrat = [];
end
%
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)
if (nverbose>=1) && shot<100000
warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check')
disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
end
return
end
systems = cellstr(systems);
if numel(systems) == 1 && strcmpi(systems{1},'all')
if shot > 24381 && shot < 49896, systems = {'main','edge'};
else, systems = {'main'};
end
end
nsys = numel(systems);
data_tmp = repmat(struct,nsys,1);
status = false(nsys,1);
for is = 1:nsys
switch lower(systems{is})
case 'main'
edge_str = '';
case 'edge'
edge_str = '.edge';
otherwise
warning('Unrecognised thomson scattering system name: %s',systems{is});
disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
return
end
nodenameeff = ['\results::thomson' edge_str ':' data_request_eff(1:2)];
tracetdi=tdi(nodenameeff);
if isempty(tracetdi.data) || ischar(tracetdi.data) || isempty(tracetdi.dim)
data_tmp(is).data = [];
data_tmp(is).error_bar = [];
data_tmp(is).data_fullpath = '';
data_tmp(is).dim = {[],[]};
data_tmp(is).system = {};
continue
end
data_tmp(is).data=tracetdi.data.'; % Thomson data as (t,z)
tracestd=tdi(['\results::thomson' edge_str ':' data_request_eff(1:2) ':error_bar']);
data_tmp(is).error_bar=tracestd.data';
data_tmp(is).data_fullpath=nodenameeff;
data_tmp(is).r_thomson = mdsdata(['\diagz::thomson_set_up' edge_str ':radial_pos']);
z=mdsdata(['\diagz::thomson_set_up' edge_str ':vertical_pos']);
data_tmp(is).dim={z,time};
data_tmp(is).units=tracetdi.units;
data_tmp(is).system=repmat({lower(systems{is})},numel(z),1);
status(is) = true;
end
% One system with good data
status_any = any(status);
isys_ref = find(status,1,'first');
gdat_data.data = vertcat(data_tmp.data);
gdat_data.error_bar = vertcat(data_tmp.error_bar);
if status_any,gdat_data.data_fullpath = [strjoin({data_tmp(status).data_fullpath},'; '),'; error_bar'];end
z = arrayfun(@(x) x.dim{1}, data_tmp, 'UniformOutput', false);
gdat_data.dim = {vertcat(z{:}),time};
gdat_data.dimunits = {'Z [m]' ; 'time [s]'};
gdat_data.x = gdat_data.dim{1};
gdat_data.t = time;
gdat_data.system=vertcat(data_tmp.system);
if status_any,gdat_data.units = data_tmp(isys_ref).units;end
if isfield(data_tmp,'r_thomson')
r_thomson = arrayfun(@(x) x.r_thomson, data_tmp, 'UniformOutput', false);
gdat_data.r_thomson = vertcat(r_thomson{:});
else
gdat_data.r_thomson = [];
end
% add fir if ne requested
if strcmp(data_request_eff(1:2),'ne')
tracefirrat_data = get_fir_thom_rat_data(shot,time,nverbose,data_option_for_firrat);
gdat_data.firrat=tracefirrat_data;
gdat_data.data_raw = gdat_data.data;
gdat_data.error_bar_raw = gdat_data.error_bar;
if status_any
gdat_data.data = gdat_data.data_raw * diag(tracefirrat_data);
gdat_data.error_bar = gdat_data.error_bar_raw * diag(tracefirrat_data);
end
gdat_data.data_fullpath=[gdat_data.data_fullpath '; fir_thom_rat ; _raw without firrat'];
if ~any(isfinite(tracefirrat_data))
gdat_data.data_fullpath=[gdat_data.data_fullpath ' WARNING use ne_raw since firrat has NaNs thus ne=NaNs; fir_thom_rat ; _raw without firrat'];
disp('***********************************************************************')
disp('WARNING: ne from Thomson has fir_thom_rat with Nans so only ne_raw can be used');
disp('***********************************************************************')
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [psiscatvol,psi_max] = get_thomson_psiscatvol(shot,systems,substr_liuqe,nverbose)
% Subfunction to join multiple systems if needed
%
systems = cellstr(systems);
if numel(systems) == 1 && strcmpi(systems{1},'all')
if shot > 24381 && shot < 49896, systems = {'main','edge'};
else, systems = {'main'};
end
end
nsys = numel(systems);
[psi_max_,psiscatvol_] = deal(cell(nsys,1));
status = false(nsys,1);
validdata = @(x) ~isempty(x.data) & ~ischar(x.data) & ~isempty(x.dim);
for is = 1:nsys
switch lower(systems{is})
case 'main'
edge_str = '';
case 'edge'
edge_str = '.edge';
otherwise
warning('Unrecognised thomson scattering system name: %s',systems{is});
disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
return
end
psi_max_{is} =gdat_tcv([],['\results::thomson',edge_str,':psi_max' substr_liuqe],'nverbose',nverbose);
psiscatvol_{is}=gdat_tcv([],['\results::thomson',edge_str,':psiscatvol' substr_liuqe],'nverbose',nverbose);
status(is) = validdata(psi_max_{is}) & validdata(psiscatvol_{is});
end
% Systems with good data
nsysv = sum(status);
if nsysv == 0
psi_max = psi_max_{1};
psiscatvol = psiscatvol_{1};
else
% Select and concatenate
psi_max_v = [psi_max_{status}];
psiscatvol_v = [psiscatvol_{status}];
psi_max = psi_max_v(1);
psi_max.data = horzcat(psi_max_v.data);
psi_max.data_fullpath = strjoin({psi_max_v.data_fullpath},'; ');
psi_max.label = strjoin({psi_max_v.data_fullpath},'; ');
psiscatvol = psiscatvol_v(1);
psiscatvol.data = horzcat(psiscatvol_v.data);
z = arrayfun(@(x) x.dim{2}, psiscatvol_v, 'UniformOutput', false);
psiscatvol.dim{2} = vertcat(z{:});
psiscatvol.t = psiscatvol.dim{2}; %TODO: Time and space are inverted for now
psiscatvol.data_fullpath = strjoin({psiscatvol_v.data_fullpath},'; ');
psiscatvol.label = strjoin({psiscatvol_v.data_fullpath},'; ');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, psitbx_str, gdat_params )
% Input arguments:
% gdat_data: current structure containing TS data
% zshift: vertical shift of the equilibrium, can be a vector for time-dependent data
% psitbx_str: source for the equilibrium map, psitbx-style
% gdat_params: parameters for the gdat call % TODO: Why is it different from gdat_data.gdat_params?
% Use psitbx
t_th = gdat_data.t;
th_z = gdat_data.x;
th_r = 0.9*ones(size(th_z));
ntt = numel(t_th);
nch = numel(th_z); % number of chords
[t_psi,status] = mdsvalue('dim_of(tcv_eq("i_p",$))',psitbx_str); % TODO: replace with time_psi when tcv_eq supports it
if ~rem(status,2)
warning('problem with getting equilibrium time base in gdat_tcv')
display(t_psi);
return
end
zshifteff=zshift;
% set zshifteff time array same as t_psi
switch numel(zshifteff)
case 1
zshifteff=zshifteff * ones(size(t_th));
case numel(t_psi)
zshifteff=interp1(t_psi,zshifteff,t_th);
case numel(t_th)
% ok
otherwise
if (gdat_params.nverbose>=1)
disp(' bad time dimension for zshift')
disp(['it should be 1 or ' num2str(numel(t_th)) ' or ' num2str(numel(t_psi))])
return
end
end
% Get LIUQE times corresponding to TS times
t_eq = t_psi(iround_os(t_psi,t_th));
% Get PSI map
psi = psitbxtcv(gdat_data.shot,t_eq,'*0',psitbx_str);
% PSITBXTCV will have removed duplicate times
i_psi = iround(psi.t,t_eq); % Mapping from psi times to TS times
% Define grid points where to evaluate PSI
if ~any(diff(zshifteff))
% Constant vector
pgrid = psitbxgrid('cylindrical','points',{th_r,th_z - zshifteff(1),0});
else
% We need as many time points in psi as TS times
psi = subs_time(psi,i_psi); % Now has ntt times
i_psi = 1:ntt; % Trivial mapping from psi times to TS times
% Time-varying vector
th_z_eff = repmat(th_z(:),1,ntt) - repmat(reshape(zshifteff,1,ntt),nch,1);
th_r_eff = repmat(th_r(:),1,ntt);
pgrid = psitbxgrid('cylindrical','time-points',{th_r_eff,th_z_eff,0});
end
% Compute psi at TS positions
psi_ts = psitbxf2f(psi,pgrid);
psiscatvol.data = squeeze(psi_ts.x(:,i_psi));
psiscatvol.dim{1} = gdat_data.x;
psiscatvol.dim{2} = t_th;
% NOTE: we should probably not include time points where equilibrium time is far from TS time.
% Compute psi_axis at TS times
% Do not use psi_axis node because:
% - For legacy LIUQE, psi_axis is not computed by psitbx, so we could get a
% different answer and cause complex numbers computing rho_psi
% - For MATLAB LIUQE, psi_axis is only the result of asxymex whereas the
% psitbx adds some Newton iterations so again complex numbers are possible
psi_norm = psitbxp2p(psi,'01');
psi_max.data = psi_norm.psimag(i_psi);
psi_max.dim = {t_th};