function [gdat_data,gdat_params,error_status,varargout] = gdat_aug(shot,data_request,varargin)
%
% function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
%
% Aim: get data from a given machine using full path or keywords.
% data_request are and should be case independent (transformed in lower case in the function and outputs)
%
% If no inputs are provided, return the list of available pre-defined data_request in gdat_data and default parameters gdat_params
%
% Inputs:
%
% no inputs: return default parameters in a structure form in gdat_params
% shot: shot number
% data_request: keyword (like 'ip') or trace name or structure containing all parameters but shot number
% varargin{i},varargin{i+1},i=1:nargin-2: additional optional parameters given in pairs: param_name, param_value
% The optional parameters list might depend on the data_request
% examples of optional parameters:
% 'plot',1 (plot is set by default to 0)
% 'machine','AUG' (the default machine is the local machine)
%
%
% Outputs:
%
% gdat_data: structure containing the data, the time trace if known and other useful information
% gdat_data.t : time trace
% gdat_data.data: requested data values
% gdat_data.dim : values of the various coordinates related to the dimensions of .data(:,:,...)
% note that one of the dim is the time, replicated in .t for clarity
% gdat_data.dimunits : units of the various dimensions, 'dimensionless' if dimensionless
% gdat_data.error_bar : if provided
% gdat_data.gdat_call : list of parameters provided in the gdat call (so can be reproduced)
% gdat_data.shot: shot number
% gdat_data.machine: machine providing the data
% gdat_data.gdat_request: keyword for gdat if relevant
% gdat_data.data_fullpath: full path to the data node if known and relevant, or relevant expression called if relevant
% gdat_data.gdat_params: copy gdat_params for completeness (gdat_params contains a .help structure detailing each parameter)
% gdat_data.xxx: any other relevant information
%
% gdat_params contains the options relevant for the called data_request. It also contains a help structure for each option
% eg.: param1 in gdat_params.param1 and help in gdat_params.help.param1
%
% Examples:
% (should add working examples for various machines (provides working shot numbers for each machine...))
%
% [a1,a2]=gdat;
% a2.data_request = 'Ip';
% a3=gdat(32827,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(32827,'ip'); % standard call
% a6=gdat(32827,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (note all lowercase in output)
% eqd = gdat(33134,'eqdsk','equil','IDE','time',2,'extra_arg_sf2sig','''-ed'',2');
%
% to get parameters, similar special signals use: 'XXX:name' with XXX in {'param:', 'param-set:', 'area-base:', 'time-base:'}
% gyrofreq=gdat(30594,{'ECS','P_sy2_g1','AUGD','param:gyr_freq'});
% or
% gyrofreq=gdat(30594,{'ECS','P_sy2_g1'},'special_signal',,'param:gyr_freq');
%
% ip=gdat(shot,'ip'); % will show in ip.gdat_fullpath the effective signal used. In this case it is similar to ask:
% ip=gdat(shot,{'FPC','IpiFP','AUGD'}); % or ip=gdat(shot,{'FPC','IpiFP'});
%
% so any signal can be loaded by giving the diagnostic, signal name and experiment, etc
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Remote data access for AUG
% You need to make a "tunnel" in one unix session from the remote node using for example:
% to create a tunnel to port 8001 to mds server lxmdsplus.aug.ipp.mpg.de
% ssh -l ipp_username -L 8001:lxmdsplus.aug.ipp.mpg.de:8000 gate1.aug.ipp.mpg.de
%
% Then in another unix session on the remote node, in matlab and with the appropriate mds path do:
% >> mdsconnect('localhost:8001')
% >> mdsvalue('1+2') % should return 3 if correctly connected
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prepare some variables, etc
varargout{1}=cell(1,1);
error_status=1;
% construct main default parameters structure
gdat_params.data_request = '';
default_machine = 'aug';
gdat_params.machine=default_machine;
gdat_params.doplot = 0;
gdat_params.exp_name = 'AUGD';
gdat_params.nverbose = 1;
% construct list of keywords from global set of keywords and specific AUG 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 AUG specific to all:
if ~isempty(data_request_names.aug)
aug_names = fieldnames(data_request_names.aug);
for i=1:length(aug_names)
data_request_names.all.(aug_names{i}) = data_request_names.aug.(aug_names{i});
end
end
data_request_names_all = 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 = aug_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_params = gdat_data.gdat_params;
% no inputs
if nargin==0
% return defaults and list of keywords
return
end
do_mdsopen_mdsclose = 1;
% treat 1st arg
if nargin>=1
if isempty(shot)
% means mdsopen(shot) already performed
shot=-999; % means not defined
do_mdsopen_mdsclose = 0;
elseif isnumeric(shot)
gdat_data.shot = shot;
elseif ischar(shot)
ivarargin_first_char = 1;
else
if gdat_params.nverbose>=1; warning('type of 1st argument unexpected, should be numeric or char'); end
error_status=2;
return
end
if nargin==1
% Only shot number given. If there is a default data_request set it and continue, otherwise return
return
end
end
% 2nd input argument if not part of pairs
if nargin>=2 && ivarargin_first_char~=1
if isempty(data_request)
return
end
% 2nd arg can be a structure with all options except shot_number, or a string for the pathname or keyword, or the start of pairs string/value for the parameters
if isstruct(data_request)
if ~isfield(data_request,'data_request')
if gdat_params.nverbose>=1; warning('expects field data_request in input parameters structure'); end
error_status=3;
return
end
%data_request.data_request = lower(data_request.data_request);
data_request_eff = data_request.data_request;
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:
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})
varargin_eff{i} = lower(varargin_eff{i});
% enforce lower case for any character driven input (only for 1st part of pair...)
% $$$ if ischar(varargin_eff{i+1})
gdat_params.(varargin_eff{i}) = varargin_eff{i+1};
% $$$ else
% $$$ gdat_params.(lower(varargin_eff{i})) = varargin_eff{i+1};
% $$$ end
else
if gdat_params.nverbose>=1; warning(['input argument nb: ' num2str(i) ' is incorrect, expects a character string']); end
error_status=401;
return
end
end
else
if gdat_params.nverbose>=1; warning('number of input arguments incorrect, cannot make pairs of parameters'); end
error_status=402;
return
end
end
data_request_eff = gdat_params.data_request; % in case was defined in pairs
% defaults for all
% default equil:
if ~isfield(gdat_params,'equil'); gdat_params.equil = 'EQI'; end
if isfield(gdat_params,'source') && (any(strcmp(lower(gdat_params.source),'ide')) || any(strcmp(lower(gdat_params.source),'idg')) ...
|| any(strcmp(lower(gdat_params.source),'ida')))
gdat_params.equil = 'IDE';
end
% $$$ if isfield(gdat_params,'source') && (any(strcmp(lower(gdat_params.source),'tre')) || any(strcmp(lower(gdat_params.source),'tra')))
% $$$ gdat_params.equil = 'TRE';
% $$$ end
extra_arg_sf2sig = '[]';
if isfield(gdat_params,'extra_arg_sf2sig') && ~isempty(gdat_params.extra_arg_sf2sig)
extra_arg_sf2sig = gdat_params.extra_arg_sf2sig;
end
gdat_params.extra_arg_sf2sig = extra_arg_sf2sig;
%
special_signal = '';
if isfield(gdat_params,'special_signal') && ~isempty(gdat_params.special_signal)
special_signal = gdat_params.special_signal;
end
gdat_params.special_signal = special_signal;
% if it is a request_keyword can obtain description:
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specifications on how to get the data provided in aug_requests_mapping
mapping_for_aug = aug_requests_mapping(data_request_eff);
gdat_data.label = mapping_for_aug.label;
% fill again at end to have full case, but here to have present status in case of early return
gdat_data.gdat_params.help = aug_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_data.mapping_for.aug = mapping_for_aug;
gdat_params = gdat_data.gdat_params;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1st treat the simplest method: "signal"
if strcmp(mapping_for_aug.method,'signal')
exp_location = gdat_data.gdat_params.exp_name;
if ~iscell(mapping_for_aug.expression)
error_status = 1010;
disp(['expects a cell array with at least 2 cells in expression, mapping_for_aug.expression = ' mapping_for_aug.expression]);
if isfield(mapping_for_aug,'not_found') && mapping_for_aug.not_found
disp(['data request not found in aug_requests_mapping, see .data_request_names_all to see all options'])
gdat_data.data_request_names_all = data_request_names_all;
end
return
elseif length(mapping_for_aug.expression)>=3 && ~isempty(mapping_for_aug.expression{3})
exp_location = mapping_for_aug.expression{3};
elseif length(mapping_for_aug.expression)>=2
mapping_for_aug.expression{3} = exp_location;
else
error_status = 101;
disp(['expects at least 2 cells in expression, mapping_for_aug.expression = ' mapping_for_aug.expression]);
if isfield(mapping_for_aug,'not_found') && mapping_for_aug.not_found
disp(['data request not found in aug_requests_mapping, see .data_request_names_all to see all options'])
gdat_data.data_request_names_all = data_request_names_all;
end
return
end
% allow changing main source with simple signals, 'source' parameter relates to mapping_for_aug.expression{1}
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || ~ischar(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = mapping_for_aug.expression{1};
else
% special case for EQI,EQH,IDE instead of IDG when FPG is present value
if strcmp(lower(gdat_data.gdat_params.source),'ide') && strcmp(lower(mapping_for_aug.expression{1}),'fpg')
gdat_data.gdat_params.source = 'IDG';
elseif strcmp(lower(gdat_data.gdat_params.source),'eqi') && strcmp(lower(mapping_for_aug.expression{1}),'fpg')
gdat_data.gdat_params.source = 'GQI';
elseif strcmp(lower(gdat_data.gdat_params.source),'eqh') && strcmp(lower(mapping_for_aug.expression{1}),'fpg')
gdat_data.gdat_params.source = 'GQH';
end
mapping_for_aug.expression{1} = gdat_data.gdat_params.source;
end
gdat_data.gdat_params.exp_name = exp_location;
% time interval
time_interval = [];
% extra args for rdaAUG_eff
if length(mapping_for_aug.expression)>=4 && ~isempty(mapping_for_aug.expression{4})
gdat_data.gdat_params.special_signal = mapping_for_aug.expression{4};
end
[aatmp,error_status]=rdaAUG_eff(shot,mapping_for_aug.expression{1},mapping_for_aug.expression{2},exp_location, ...
time_interval,gdat_data.gdat_params.extra_arg_sf2sig,gdat_data.gdat_params.special_signal);
if error_status~=0
if gdat_params.nverbose>=3; disp(['error after rdaAUG in signal with data_request_eff= ' data_request_eff]); end
return
end
if isfield(aatmp,'data'); gdat_data.data = aatmp.data; end
if isfield(aatmp,'t'); gdat_data.t = aatmp.t; end
if isfield(aatmp,'x'); gdat_data.x = aatmp.x; end
if isfield(aatmp,'dimunits'); gdat_data.dimunits = aatmp.dimunits; end
gdat_data.data_fullpath=mapping_for_aug.expression;
if isempty(aatmp.data) || (isempty(gdat_data.t) && isempty(gdat_data.x))
return
end
if mapping_for_aug.timedim<=0
% need to guess timedim
if prod(size(aatmp.data)) == length(aatmp.data)
mapping_for_aug.timedim = 1;
elseif length(size(aatmp.data))==2
% 2 true dimensions
if length(aatmp.t) == size(aatmp.data,1)
mapping_for_aug.timedim = 1;
else
mapping_for_aug.timedim = 2;
end
else
% more than 2 dimensions, find the one with same length to define timedim
mapping_for_aug.timedim = find(size(aatmp.data)==length(aatmp.t));
if length(mapping_for_aug.timedim); mapping_for_aug.timedim = mapping_for_aug.timedim(1); end
end
mapping_for_aug.gdat_timedim = mapping_for_aug.timedim;
end
if length(size(aatmp.data))==1 | (length(size(aatmp.data))==2 & size(aatmp.data,2)==1)
gdat_data.dim=[{aatmp.t}];
gdat_data.dimunits={'time [s]'};
elseif length(size(aatmp.data))==2
gdat_data.dim{mapping_for_aug.timedim} = gdat_data.t;
gdat_data.dimunits{mapping_for_aug.timedim} = 'time [s]';
gdat_data.dim{setdiff([1:2],mapping_for_aug.timedim)} = gdat_data.x;
else
gdat_data.dim{mapping_for_aug.timedim} = gdat_data.t;
gdat_data.dimunits{mapping_for_aug.timedim} = 'time [s]';
nbdims = length(size(gdat_data.data));
for i=1:nbdims
if i~=mapping_for_aug.timedim
ieff=i;
if i>mapping_for_aug.timedim; ieff=i-1; end
gdat_data.dim{i} = gdat_data.x{ieff};
end
end
end
nbdims = length(gdat_data.dim);
gdat_data.units = aatmp.units;
if mapping_for_aug.gdat_timedim>0 && mapping_for_aug.gdat_timedim ~= mapping_for_aug.timedim
% shift timedim to gdat_timedim data(i,j,...itime,k,...) -> data(i,inewtime,j,...,k,...)
% note that this means that gdat_data.x and gdat_data.t are same and correct,
% only .data, .dim and .dimunits need to be changed
iprev=[1:nbdims];
ij=find(dim_nontim>mapping_for_aug.gdat_timedim-1);
inew=[1:mapping_for_aug.gdat_timedim-1 mapping_for_aug.timedim dim_nontim(ij)];
data_sizes = size(aatmp.data);
gdat_data.data = NaN*ones(data_sizes(inew));
abcol=ones(1,nbdims)*double(':'); abcomma=ones(1,nbdims)*double(',');
dimstr_prev=['(' repmat(':,',1,mapping_for_aug.timedim-1) 'it,' ...
repmat(':,',1,nbdims-mapping_for_aug.timedim-1) ':)'];
dimstr_new=['(' repmat(':,',1,mapping_for_aug.gdat_timedim-1) 'it,' ...
repmat(':,',1,nbdims-mapping_for_aug.gdat_timedim-1) ':)'];
% eval gdat_data.data(;,:,...,it,...) = aatmp.data(:,:,:,it,...);
for it=1:size(aatmp.data,mapping_for_aug.timedim)
shift_eval = ['gdat_data.data' dimstr_new ' = aatmp.data' dimstr_prev ';'];
eval(shift_eval);
end
gdat_data.dim = aatmp.dim(inew);
% gdat_data.dimunits = aatmp.dimunits(inew);
else
mapping_for_aug.gdat_timedim = mapping_for_aug.timedim;
end
% end of method "signal"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_aug.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_aug.expression ';'];
gdata_data_orig = gdat_data;
fields_to_copy = {'gdat_params', 'gdat_request', 'label', 'data_fullpath', 'request_description', 'mapping_for'};
eval([mapping_for_aug.expression ';']);
for i=1:length(fields_to_copy)
if isfield(gdata_data_orig,fields_to_copy{i})
gdat_tmp.(fields_to_copy{i}) = gdata_data_orig.(fields_to_copy{i});
end
end
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_aug.expression]); end
error_status=801;
return
end
tmp_fieldnames = fieldnames(gdat_tmp);
if sum(strcmp(tmp_fieldnames,'data'))==0 % note: cannot do isfield since gdat_tmp might be an object
if (gdat_params.nverbose>=1); warning(['expression does not return a child name ''data'' for ' data_request_eff]); end
end
for i=1:length(tmp_fieldnames)
gdat_data.(tmp_fieldnames{i}) = gdat_tmp.(tmp_fieldnames{i});
end
% add .t and .x in case only dim is provided
% do not allow shifting of timedim since should be treated in the relevant expression
ijdim=find(strcmp(tmp_fieldnames,'dim')==1);
if ~isempty(ijdim)
nbdims = length(gdat_data.dim);
if mapping_for_aug.timedim==-1;
mapping_for_aug.timedim = nbdims;
if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_aug.timedim = nbdims-1; end
end
dim_nontim = setdiff([1:nbdims],mapping_for_aug.timedim);
ijt=find(strcmp(tmp_fieldnames,'t')==1);
if isempty(ijt)
gdat_data.t = gdat_data.dim{mapping_for_aug.timedim};
end
ijx=find(strcmp(tmp_fieldnames,'x')==1);
if isempty(ijx)
if ~isempty(dim_nontim)
% since most cases have at most 2d, copy as array if data is 2D and as cell if 3D or more
if length(dim_nontim)==1
gdat_data.x = gdat_data.dim{dim_nontim(1)};
else
gdat_data.x = gdat_data.dim(dim_nontim);
end
end
end
gdat_data.data_fullpath=mapping_for_aug.expression;
end
% end of method "expression"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(mapping_for_aug.method,'switchcase')
switch data_request_eff % not lower(...) since data_request_eff should be lower case already at this stage
case {'cxrs', 'cxrs_rho'}
% if 'fit' option is added: 'fit',1, then the fitted profiles are returned which means "_rho" is effectively called
if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit) || ~isnumeric(gdat_data.gdat_params.fit)
if strcmp(data_request_eff,'cxrs')
gdat_data.gdat_params.fit = 0;
else
gdat_data.gdat_params.fit = 1; % add fit as default if cxrs_rho requested (since equil takes most time)
end
elseif gdat_data.gdat_params.fit==1
data_request_eff = 'cxrs_rho';
end
exp_name_eff = gdat_data.gdat_params.exp_name;
sources_available = {'CEZ', 'CMZ', 'CUZ', 'COZ'};
r_node_available = {'R_time', 'R', 'R_time', 'R_time'};
z_node_available = {'z_time', 'z', 'z_time', 'z_time'};
% scroll through list up to first available when no sources provided
sources_available_for_scroll = {'CEZ', 'CUZ', 'COZ'};
scroll_through_list = 0;
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = sources_available_for_scroll(1);
scroll_through_list = 1; % means scroll through sources until one is available
elseif ischar(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = {gdat_data.gdat_params.source};
end
if length(gdat_data.gdat_params.source)==1
if strcmp(upper(gdat_data.gdat_params.source{1}),'CEZ')
gdat_data.gdat_params.source = {'CEZ'};
elseif strcmp(upper(gdat_data.gdat_params.source{1}),'CMZ')
gdat_data.gdat_params.source = {'CMZ'};
elseif strcmp(upper(gdat_data.gdat_params.source{1}),'CUZ')
gdat_data.gdat_params.source = {'CUZ'};
elseif strcmp(upper(gdat_data.gdat_params.source{1}),'COZ')
gdat_data.gdat_params.source = {'COZ'};
elseif strcmp(upper(gdat_data.gdat_params.source{1}),'ALL')
gdat_data.gdat_params.source = sources_available;
scroll_through_list = 2; % means scroll through all sources and load all sources
else
warning(['source = ' gdat_data.gdat_params.source ' not expected with data_request= ' data_request_eff]);
return
end
else
sources_in = intersect(sources_available,upper(gdat_data.gdat_params.source));
if length(sources_in) ~= length(gdat_data.gdat_params.source)
disp('following sources not yet available, check with O. Sauter if need be')
setdiff(upper(gdat_data.gdat_params.source),sources_available)
end
gdat_data.gdat_params.source = sources_in;
end
extra_arg_sf2sig_eff_string = '';
if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
extra_arg_sf2sig_eff_string = [',' gdat_data.gdat_params.extra_arg_sf2sig];
end
% set starting source
i_count = 1;
diag_name = gdat_data.gdat_params.source{i_count};
sources_tried{i_count} = diag_name;
iload = 1;
iequil = 0;
while iload==1
ishotfile_ok = 1;
i_source = strmatch(diag_name,sources_available);
r_node = r_node_available{i_source};
z_node = z_node_available{i_source};
% R, Z positions of measurements
try
% eval(['[r_time]=sf2ab(diag_name,shot,r_node,''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
% both for CEZ and CMZ, and.. Ti:1 is ok, otherwise introduce string above
[r_time,err]=rdaAUG_eff(shot,diag_name,r_node,exp_name_eff,[],extra_arg_sf2sig_eff_string,['area-base:' r_node ':1']);
catch ME_R_time
% assume no shotfile
disp(getReport(ME_R_time))
ishotfile_ok = 0;
end
if ishotfile_ok == 1
gdat_data.r = r_time.data;
inotok=find(gdat_data.r<=0);
gdat_data.r(inotok) = NaN;
try
% eval(['[z_time]=sf2ab(diag_name,shot,z_node,''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
[z_time,err]=rdaAUG_eff(shot,diag_name,z_node,exp_name_eff,[],extra_arg_sf2sig_eff_string,['area-base:' z_node ':1']);
gdat_data.z = z_time.data;
inotok=find(gdat_data.z<=0);
gdat_data.z(inotok) = NaN;
catch ME_R_time
disp(getReport(ME_R_time))
ishotfile_ok = 0;
end
else
gdat_data.r = [];
gdat_data.z = [];
end
if ishotfile_ok == 1
try
% eval(['[time]=sf2tb(diag_name,shot,''time'',''-exp'',exp_name_eff' extra_arg_sf2sig_eff_string ');']);
[time,err]=rdaAUG_eff(shot,diag_name,'time',exp_name_eff,[],extra_arg_sf2sig_eff_string,'time-base:time:0');
gdat_data.t = time.data;
catch ME_R_time
disp(getReport(ME_R_time))
ishotfile_ok = 0;
end
else
gdat_data.t = [];
end
gdat_data.dim{1} = {gdat_data.r , gdat_data.z};
gdat_data.dimunits{1} = 'R, Z [m]';
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits{2} = 't [s]';
gdat_data.x = gdat_data.dim{1};
% vrot
if ishotfile_ok == 1
try
[a,error_status]=rdaAUG_eff(shot,diag_name,'vrot',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
if isempty(a.data) || isempty(a.t) || error_status>0
if gdat_params.nverbose>=3;
a
disp(['with data_request= ' data_request_eff])
end
end
[aerr,e]=rdaAUG_eff(shot,diag_name,'err_vrot',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
catch ME_vrot
disp(getReport(ME_vrot))
ishotfile_ok = 0;
end
gdat_data.vrot.data = a.data;
gdat_data.vrot.error_bar = aerr.data;
else
gdat_data.vrot.data = [];
gdat_data.vrot.error_bar = [];
end
a.name = 'vrot';
if any(strcmp(fieldnames(a),'units')); gdat_data.vrot.units=a.units; end
if any(strcmp(fieldnames(a),'name')); gdat_data.vrot.name=a.name; end
gdat_data.vrot.label = 'vrot_tor';
% Ti
% [a,e]=rdaAUG_eff(shot,diag_name,'Ti',exp_name_eff);
% [aerr,e]=rdaAUG_eff(shot,diag_name,'err_Ti',exp_name_eff);
if ishotfile_ok == 1
try
[a,e]=rdaAUG_eff(shot,diag_name,'Ti_c',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
[aerr,e]=rdaAUG_eff(shot,diag_name,'err_Ti_c',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
catch ME_ti
disp(getReport(ME_ti))
ishotfile_ok = 0;
end
gdat_data.ti.data = a.data;
gdat_data.ti.error_bar = aerr.data;
else
gdat_data.ti.data = [];
gdat_data.ti.error_bar = [];
end
a.name = 'Ti_c';
if any(strcmp(fieldnames(a),'units')); gdat_data.ti.units=a.units; end
if any(strcmp(fieldnames(gdat_data.ti),'units')); gdat_data.units=gdat_data.ti.units; end
if any(strcmp(fieldnames(a),'name')); gdat_data.ti.name=a.name; end
gdat_data.ti.label = 'Ti_c';
% main node ti a this stage
gdat_data.data = gdat_data.ti.data;
gdat_data.label = [gdat_data.label '/Ti'];
gdat_data.error_bar = gdat_data.ti.error_bar;
gdat_data.data_fullpath=[exp_name_eff '/' diag_name '/' a.name ';vrot, Ti in data, {r;z} in .x'];
%
if strcmp(data_request_eff,'cxrs_rho')
% defaults
if iequil == 0
gdat_data.equil.data = [];
gdat_data.psi = [];
gdat_data.rhopolnorm = [];
gdat_data.rhotornorm = [];
gdat_data.rhovolnorm = [];
% defaults for fits, so user always gets std structure
gdat_data.fit.rhotornorm = []; % same for both ti and vrot
gdat_data.fit.rhopolnorm = [];
gdat_data.fit.t = [];
gdat_data.fit.ti.data = [];
gdat_data.fit.ti.drhotornorm = [];
gdat_data.fit.vrot.data = [];
gdat_data.fit.vrot.drhotornorm = [];
gdat_data.fit.raw.rhotornorm = [];
gdat_data.fit.raw.ti.data = [];
gdat_data.fit.raw.vrot.data = [];
fit_tension_default = -1;
if isfield(gdat_data.gdat_params,'fit_tension')
fit_tension = gdat_data.gdat_params.fit_tension;
else
fit_tension = fit_tension_default;
end
if ~isstruct(fit_tension)
fit_tension_eff.ti = fit_tension;
fit_tension_eff.vrot = fit_tension;
fit_tension = fit_tension_eff;
else
if ~isfield(fit_tension,'ti'); fit_tension.ti = fit_tension_default; end
if ~isfield(fit_tension,'vrot '); fit_tension.vrot = fit_tension_default; end
end
gdat_data.gdat_params.fit_tension = fit_tension;
if isfield(gdat_data.gdat_params,'fit_nb_rho_points')
fit_nb_rho_points = gdat_data.gdat_params.fit_nb_rho_points;
else
fit_nb_rho_points = 201;
end
gdat_data.gdat_params.fit_nb_rho_points = fit_nb_rho_points;
end
if ishotfile_ok == 1 && iequil == 0
params_equil = gdat_data.gdat_params;
params_equil.data_request = 'equil';
[equil,params_equil,error_status] = gdat_aug(shot,params_equil);
if error_status>0
if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end
return
end
iequil = 1;
gdat_data.gdat_params.equil = params_equil.equil;
gdat_data.equil = equil;
inb_chord_cxrs=size(gdat_data.data,1);
inb_time_cxrs=size(gdat_data.data,2);
psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhopolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[min(gdat_data.t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.t(end)+0.1)];
iok=find(~isnan(gdat_data.r(:,1)));
for itequil=1:length(time_equil)-1
rr=equil.Rmesh(:,itequil);
zz=equil.Zmesh(:,itequil);
psirz_in = equil.psi2D(:,:,itequil);
it_cxrs_inequil = find(gdat_data.t>time_equil(itequil) & gdat_data.t<=time_equil(itequil+1));
if ~isempty(it_cxrs_inequil)
if ~strcmp(upper(gdat_data.gdat_params.source),'CMZ')
rout=gdat_data.r(iok);
zout=gdat_data.z(iok);
else
rout=gdat_data.r(iok,it_cxrs_inequil);
zout=gdat_data.z(iok,it_cxrs_inequil);
end
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
if ~strcmp(upper(gdat_data.gdat_params.source),'CMZ')
psi_out(iok,it_cxrs_inequil) = repmat(psi_at_routzout,[1,length(it_cxrs_inequil)]);
else
psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil));
end
rhopolnorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
for it_cx=1:length(it_cxrs_inequil)
rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
end
gdat_data.psi = psi_out;
gdat_data.rhopolnorm = rhopolnorm_out;
gdat_data.rhotornorm = rhotornorm_out;
gdat_data.rhovolnorm = rhovolnorm_out;
aaa.shot = gdat_data.shot;
aaa.data = gdat_data.data;
aaa.x = gdat_data.rhopolnorm;
aaa.t = gdat_data.t;
aaa.gdat_params = gdat_data.gdat_params;
aaa = get_grids_1d(aaa,2,1);
gdat_data.grids_1d = aaa.grids_1d;
%
if gdat_data.gdat_params.fit==1
% add fits
gdat_data.fit.raw.rhotornorm = NaN*ones(size(gdat_data.ti.data));
gdat_data.fit.raw.ti.data = NaN*ones(size(gdat_data.ti.data));
gdat_data.fit.raw.vrot.data = NaN*ones(size(gdat_data.vrot.data));
rhotornormfit = linspace(0,1,fit_nb_rho_points)';
gdat_data.fit.rhotornorm = rhotornormfit;
gdat_data.fit.t = gdat_data.t;
for it=1:length(gdat_data.t)
% make rhotor->rhopol transformation for each time since equilibrium might have changed
gdat_data.fit.rhopolnorm(:,it)=interpos(gdat_data.grids_1d.rhotornorm_equil(:,it), ...
gdat_data.grids_1d.rhopolnorm_equil(:,it),rhotornormfit);
idata = find(gdat_data.ti.data(:,it)>0 & gdat_data.rhotornorm(:,it)<1.01);
if length(idata)>0
gdat_data.fit.ti.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it);
gdat_data.fit.ti.raw.data(idata,it) = gdat_data.ti.data(idata,it);
gdat_data.fit.ti.raw.error_bar(idata,it) = gdat_data.ti.error_bar(idata,it);
gdat_data.fit.vrot.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it);
gdat_data.fit.vrot.raw.data(idata,it) = gdat_data.vrot.data(idata,it);
gdat_data.fit.vrot.raw.error_bar(idata,it) = gdat_data.vrot.error_bar(idata,it);
[rhoeff,irhoeff] = sort(gdat_data.rhotornorm(idata,it));
rhoeff = [0; rhoeff];
% they are some strange low error_bars, so remove these by max(Ti/error_bar)<=10; and changing it to large error bar
tieff = gdat_data.ti.data(idata(irhoeff),it);
ti_err_eff = gdat_data.ti.error_bar(idata(irhoeff),it);
ij=find(tieff./ti_err_eff>10.);
if ~isempty(ij); ti_err_eff(ij) = tieff(ij)./0.1; end
vroteff = gdat_data.vrot.data(idata(irhoeff),it);
vrot_err_eff = gdat_data.vrot.error_bar(idata(irhoeff),it);
ij=find(vroteff./vrot_err_eff>10.);
if ~isempty(ij); vrot_err_eff(ij) = vroteff(ij)./0.1; end
%
tieff = [tieff(1); tieff];
ti_err_eff = [1e4; ti_err_eff];
vroteff = [vroteff(1); vroteff];
vrot_err_eff = [1e5; vrot_err_eff];
[gdat_data.fit.ti.data(:,it), gdat_data.fit.ti.drhotornorm(:,it)] = interpos(rhoeff,tieff,rhotornormfit,fit_tension.ti,[1 0],[0 0],ti_err_eff);
[gdat_data.fit.vrot.data(:,it), gdat_data.fit.vrot.drhotornorm(:,it)] = interpos(rhoeff,vroteff,rhotornormfit,fit_tension.vrot,[1 0],[0 0],vrot_err_eff);
end
end
end
end
end
if scroll_through_list == 0 || scroll_through_list == 2
if scroll_through_list == 2
tmp.(diag_name) = gdat_data;
end
if length(gdat_data.gdat_params.source) > i_count
i_count = i_count + 1;
diag_name = gdat_data.gdat_params.source{i_count};
else
iload = 0;
end
elseif scroll_through_list == 1
if ishotfile_ok == 1
iload = 0;
else
sources_remaining = setdiff(sources_available_for_scroll,sources_tried,'stable');
if ~isempty(sources_remaining)
i_count = i_count + 1;
diag_name = sources_remaining{1};
sources_tried{i_count} = diag_name;
else
iload = 0;
end
end
else
disp('should not arrive here, check value of scroll_through_list')
scroll_through_list
iload = 0;
end
end
if scroll_through_list == 2
tmp_field=fieldnames(tmp);
for i=1:length(tmp_field)
gdat_data.(tmp_field{i}) = tmp.(tmp_field{i});
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ece', 'eced', 'ece_rho', 'eced_rho'}
nth_points = 13;
if isfield(gdat_data.gdat_params,'nth_points') && ~isempty(gdat_data.gdat_params.nth_points)
nth_points = gdat_data.gdat_params.nth_points;
else
gdat_data.gdat_params.nth_points = nth_points;
end
channels = -1;
if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels)
channels = gdat_data.gdat_params.channels;
end
if nth_points>=10
match_rz_to_time = 1;
else
match_rz_to_time = 0;
end
if isfield(gdat_data.gdat_params,'match_rz_to_time') && ~isempty(gdat_data.gdat_params.match_rz_to_time)
match_rz_to_time = gdat_data.gdat_params.match_rz_to_time;
else
gdat_data.gdat_params.match_rz_to_time = match_rz_to_time;
end
time_interval = [];
if isfield(gdat_data.gdat_params,'time_interval') && ~isempty(gdat_data.gdat_params.time_interval)
time_interval = gdat_data.gdat_params.time_interval;
else
gdat_data.gdat_params.time_interval = time_interval;
end
%
if isfield(gdat_data.gdat_params,'diag_name') && ~isempty(gdat_data.gdat_params.diag_name)
diag_name = gdat_data.gdat_params.diag_name;
if strcmp(upper(diag_name),'RMD'); gdat_data.gdat_params.exp_name = 'ECED'; end
else
diag_name = 'CEC';
gdat_data.gdat_params.diag_name = diag_name;
end
exp_name_eff = gdat_data.gdat_params.exp_name;
if strcmp(data_request_eff,'eced')
exp_name_eff = 'ECED';
gdat_data.gdat_params.exp_name = exp_name_eff;
end
[a,e]=rdaAUG_eff(shot,diag_name,'Trad-A',exp_name_eff,time_interval,gdat_data.gdat_params.extra_arg_sf2sig);
% [a,e]=rdaAUG_eff(shot,diag_name,'Trad-A',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
inb_chord = size(a.data,1);
if channels(1)<=0
channels = [1:inb_chord];
end
gdat_data.dim{1} = channels;
gdat_data.gdat_params.channels = channels;
if nth_points>1
gdat_data.data = a.data(channels,1:nth_points:end);
gdat_data.dim{2} = a.t(1:nth_points:end);
else
gdat_data.data = a.data(channels,:);
gdat_data.dim{2} = a.t;
end
gdat_data.x = gdat_data.dim{1};
gdat_data.t = gdat_data.dim{2};
gdat_data.dimunits=[{'channels'} ; {'time [s]'}];
if any(strcmp(fieldnames(a),'units')); gdat_data.units=a.units; end
try
[aR,e]=rdaAUG_eff(shot,diag_name,'R-A',exp_name_eff,time_interval,gdat_data.gdat_params.extra_arg_sf2sig);
catch
end
try
[aZ,e]=rdaAUG_eff(shot,diag_name,'z-A',exp_name_eff,time_interval,gdat_data.gdat_params.extra_arg_sf2sig);
catch
disp(['problem with getting z-A in ' diag_name])
end
if match_rz_to_time
% interpolate R structure on ece data time array, to ease plot vs R
for i=1:length(channels)
radius.data(i,:) = interp1(aR.t,aR.data(channels(i),:),gdat_data.t);
zheight.data(i,:) = interp1(aZ.t,aZ.data(channels(i),:),gdat_data.t);
end
radiuszheight.t=gdat_data.t;
else
radius.data = aR.data(channels,:);
radiuszheight.t=aR.t;
zheight.data = aZ.data(channels,:);
end
ij=find(gdat_data.data<=0);
gdat_data.data(ij)=NaN;
gdat_data.rz.r=radius.data;
gdat_data.rz.z=zheight.data;
gdat_data.rz.t = radiuszheight.t;
gdat_data.data_fullpath = [exp_name_eff '/' diag_name '/Trad-A with R, Z in .r,.z'];
if strcmp(data_request_eff,'ece_rho') || strcmp(data_request_eff,'eced_rho')
params_equil = gdat_data.gdat_params;
params_equil.data_request = 'equil';
[equil,params_equil,error_status] = gdat_aug(shot,params_equil);
if error_status>0
if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end
return
end
gdat_data.gdat_params.equil = params_equil.equil;
gdat_data.equil = equil;
gdat_data.data_fullpath = [exp_name_eff '/' diag_name '/Trad-A with .r,.z projected on equil ' gdat_data.gdat_params.equil ' in .rhos'];
inb_chord_ece=size(gdat_data.rz.r,1);
inb_time_ece=size(gdat_data.rz.r,2);
psi_out = NaN*ones(inb_chord_ece,inb_time_ece);
rhopolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
rhotornorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
rhovolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[min(gdat_data.rz.t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.rz.t(end)+0.1)];
for itequil=1:length(time_equil)-1
rr=equil.Rmesh(:,itequil);
zz=equil.Zmesh(:,itequil);
psirz_in = equil.psi2D(:,:,itequil);
it_ece_inequil = find(gdat_data.rz.t>time_equil(itequil) & gdat_data.rz.t<=time_equil(itequil+1));
if ~isempty(it_ece_inequil)
r_it_ece_inequil = gdat_data.rz.r(:,it_ece_inequil);
iok = find(~isnan(r_it_ece_inequil) & r_it_ece_inequil>0);
if ~isempty(iok)
rout=r_it_ece_inequil(iok);
z_it_ece_inequil = gdat_data.rz.z(:,it_ece_inequil);
zout=z_it_ece_inequil(iok);
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
[ieff,jeff]=ind2sub(size(gdat_data.rz.r(:,it_ece_inequil)),iok);
for ij=1:length(iok)
psi_out(ieff(ij),it_ece_inequil(jeff(ij))) = psi_at_routzout(ij);
end
rhopolnorm_out(:,it_ece_inequil) = sqrt(abs((psi_out(:,it_ece_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))));
for it_cx=1:length(it_ece_inequil)
rhotornorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
rhovolnorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
end
end
gdat_data.rhos.psi = psi_out;
gdat_data.rhos.rhopolnorm = rhopolnorm_out;
gdat_data.rhos.rhotornorm = rhotornorm_out;
gdat_data.rhos.rhovolnorm = rhovolnorm_out;
gdat_data.rhos.t = gdat_data.rz.t;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'eqdsk'}
%
time_eqdsks=1.; % default time
if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time)
time_eqdsks = gdat_data.gdat_params.time;
else
gdat_data.gdat_params.time = time_eqdsks;
disp(['"time" is expected as an option, choose default time = ' num2str(time_eqdsks)]);
end
gdat_data.gdat_params.time = time_eqdsks;
gdat_data.t = time_eqdsks;
zshift = 0.;
if isfield(gdat_data.gdat_params,'zshift') && ~isempty(gdat_data.gdat_params.zshift)
zshift = gdat_data.gdat_params.zshift;
else
gdat_data.gdat_params.zshift = zshift;
end
gdat_data.gdat_params.zshift = zshift;
dowrite = 1;
if isfield(gdat_data.gdat_params,'write') && ~isempty(gdat_data.gdat_params.write)
dowrite = gdat_data.gdat_params.write;
else
gdat_data.gdat_params.write = dowrite;
end
gdat_data.gdat_params.write = dowrite;
params_equil = gdat_data.gdat_params;
params_equil.data_request = 'equil';
[equil,params_equil,error_status] = gdat_aug(shot,params_equil);
if error_status>0
if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end
return
end
gdat_data.gdat_params = params_equil;
gdat_data.equil = equil;
if gdat_data.gdat_params.doplot==0
fignb_handle = -1;
else
figure(9999);
fignb_handle = gcf;
end
for itime=1:length(time_eqdsks)
time_eff = time_eqdsks(itime);
% use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2
[eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(equil,time_eff,gdat_data.gdat_params.zshift,'source',gdat_data.gdat_params.equil,'extra_arg_sf2sig',gdat_data.gdat_params.extra_arg_sf2sig,'fignb',fignb_handle);
eqdskAUG.fnamefull = fullfile(['/tmp/' getenv('USER')],['EQDSK_' num2str(shot) 't' num2str(time_eff,'%.4f')]);
cocos_out = equil.cocos;
if isfield(gdat_data.gdat_params,'cocos') && ~isempty(gdat_data.gdat_params.cocos)
cocos_out = gdat_data.gdat_params.cocos;
else
gdat_data.gdat_params.cocos = cocos_out;
end
if equil.cocos ~= cocos_out
[eqdsk_cocosout, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdskAUG,[equil.cocos cocos_out]);
else
eqdsk_cocosout = eqdskAUG;
end
% for several times, use array of structure for eqdsks,
% cannot use it for psi(R,Z) in .data and .x since R, Z might be different at different times,
% so project psi(R,Z) on Rmesh, Zmesh of 1st time
if length(time_eqdsks) > 1
gdat_data.eqdsk{itime} = eqdsk_cocosout;
if gdat_data.gdat_params.write; gdat_data.eqdsk{itime} = write_eqdsk(eqdskAUG.fnamefull,eqdsk_cocosout); end
if itime==1
gdat_data.data(:,:,itime) = gdat_data.eqdsk{itime}.psi;
gdat_data.dim{1} = gdat_data.eqdsk{itime}.rmesh;
gdat_data.dim{2} = gdat_data.eqdsk{itime}.zmesh;
else
xx=repmat(reshape(gdat_data.dim{1},length(gdat_data.dim{1}),1),1,size(gdat_data.eqdsk{itime}.psi,2));
yy=repmat(reshape(gdat_data.dim{2},1,length(gdat_data.dim{2})),size(gdat_data.eqdsk{itime}.psi,1),1);
aa = interpos2Dcartesian(gdat_data.eqdsk{itime}.rmesh,gdat_data.eqdsk{itime}.zmesh ...
,gdat_data.eqdsk{itime}.psi,xx,yy,-1,-1);
gdat_data.data(:,:,itime) = aa;
end
else
gdat_data.eqdsk = eqdsk_cocosout;
if gdat_data.gdat_params.write; gdat_data.eqdsk = write_eqdsk(eqdskAUG.fnamefull,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);
gdat_data.data_fullpath=['psi(R,Z) and eqdsk from geteqdskAUG with ' gdat_data.gdat_params.equil ' ;zshift=' num2str(zshift)];
gdat_data.units = 'T m^2';
gdat_data.dimunits = {'m','m','s'};
gdat_data.request_description = ['data=psi, x=(R,Z), eqdsk contains eqdsk structure with which ' ...
'plot_eqdsk, write_eqdsk, read_eqdsk can be used'];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'equil'}
% get equil params and time array in gdat_data.t
[gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data);
if isempty(Lpf1_t)
disp('Lpf1_t is empty, probably no data, return')
return
end
% since Lpf depends on time, need to load all first and then loop over time for easier mapping
[qpsi,e]=rdaAUG_eff(shot,DIAG,'Qpsi',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
ndimrho = size(qpsi.data,2);
if ndimrho==NTIME_Lpf
% data seems to be transposed
ndimrho = size(qpsi.data,1);
itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
else
itotransposeback = 0;
end
qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback);
ijnan=find(isnan(qpsi.value));
qpsi.value(ijnan)=0;
[psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
[phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
[Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
[Area,e]=rdaAUG_eff(shot,DIAG,'Area',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
Area=adapt_rda(Area,NTIME,2*ndimrho,itotransposeback);
[Ri,e]=rdaAUG_eff(shot,DIAG,'Ri',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
Ri=adapt_rda(Ri,NTIME,M_Rmesh,itotransposeback);
[Zj,e]=rdaAUG_eff(shot,DIAG,'Zj',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
Zj=adapt_rda(Zj,NTIME,N_Zmesh,itotransposeback);
if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
[PFM_tree,e]=rdaAUG_eff(shot,DIAG,'PFM',exp_name_eff,[],['''-raw'',' gdat_data.gdat_params.extra_arg_sf2sig]); % -raw necessary for IDE
else
[PFM_tree,e]=rdaAUG_eff(shot,DIAG,'PFM',exp_name_eff,[],['''-raw''']); % -raw necessary for IDE
end
PFM_tree=adaptPFM_rda(PFM_tree,M_Rmesh,N_Zmesh,NTIME);
[Pres,e]=rdaAUG_eff(shot,DIAG,'Pres',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
Pres=adapt_rda(Pres,NTIME,2*ndimrho,itotransposeback);
[Jpol,e]=rdaAUG_eff(shot,DIAG,'Jpol',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
Jpol=adapt_rda(Jpol,NTIME,2*ndimrho,itotransposeback);
[FFP,e]=rdaAUG_eff(shot,DIAG,'FFP',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
if ~isempty(FFP.value)
FFP=adapt_rda(FFP,NTIME,ndimrho,itotransposeback);
else
FFP.value=NaN*ones(NTIME,max(Lpf1_t));
end
if strcmp(DIAG,'EQI') || strcmp(DIAG,'EQH') || strcmp(DIAG,'IDE')
[Rinv,e]=rdaAUG_eff(shot,DIAG,'Rinv',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
Rinv=adapt_rda(Rinv,NTIME,ndimrho,itotransposeback);
[R2inv,e]=rdaAUG_eff(shot,DIAG,'R2inv',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
R2inv=adapt_rda(R2inv,NTIME,ndimrho,itotransposeback);
[Bave,e]=rdaAUG_eff(shot,DIAG,'Bave',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
Bave=adapt_rda(Bave,NTIME,ndimrho,itotransposeback);
[B2ave,e]=rdaAUG_eff(shot,DIAG,'B2ave',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
B2ave=adapt_rda(B2ave,NTIME,ndimrho,itotransposeback);
if strcmp(DIAG,'IDE')
FTRA.value=[];
else
[FTRA,e]=rdaAUG_eff(shot,DIAG,'FTRA',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
FTRA=adapt_rda(FTRA,NTIME,ndimrho,itotransposeback);
end
else
Rinv.value=[]; R2inv.value=[]; Bave.value=[]; B2ave.value=[]; FTRA.value=[];
end
[LPFx,e]=rdaAUG_eff(shot,DIAG,'LPFx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
LPFx.value=LPFx.value(1:NTIME); LPFx.data=LPFx.value(1:NTIME); LPFx.t=LPFx.t(1:NTIME);
[PFxx,e]=rdaAUG_eff(shot,DIAG,'PFxx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
PFxx=adapt_rda(PFxx,NTIME,max(LPFx.value)+1,itotransposeback);
[RPFx,e]=rdaAUG_eff(shot,DIAG,'RPFx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
RPFx=adapt_rda(RPFx,NTIME,max(LPFx.value)+1,itotransposeback);
[zPFx,e]=rdaAUG_eff(shot,DIAG,'zPFx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
zPFx=adapt_rda(zPFx,NTIME,max(LPFx.value)+1,itotransposeback);
% seems "LCFS" q-value is far too large, limit to some max (when diverted)
max_qValue = 30.; % Note could just put a NaN on LCFS value since ill-defined when diverted
for it=1:NTIME
Lpf1 = Lpf1_t(it);
% Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
% change it to (radial,time) and use only Lpf+1 points up to LCFS
ijok=find(abs(qpsi.value)>0); % note: eqr fills in only odd points radially
% set NaNs to zeroes
if qpsi.value(ijok(1))<0
gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
else
gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
end
% get x values
gdat_data.psi(:,it)=psi_tree.value(it,Lpf1:-1:1)';
gdat_data.psi_axis(it)= gdat_data.psi(1,it);
gdat_data.psi_lcfs(it)= gdat_data.psi(end,it);
gdat_data.rhopolnorm(:,it) = sqrt(abs((gdat_data.psi(:,it)-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it))));
if strcmp(DIAG,'EQR');
% q value has only a few values and from center to edge, assume they are from central rhopol values on
% But they are every other point starting from 3rd
ijk=find(gdat_data.qvalue(:,it)~=0);
if length(ijk)>2
% now shots have non-zero axis values in eqr
rhoeff=gdat_data.rhopolnorm(ijk,it);
qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above
if ijk(1)>1
rhoeff = [0.; rhoeff];
qeff = [qeff(1) ;qeff];
end
ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it)));
qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]);
else
qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
end
gdat_data.qvalue(:,it) = qfit;
end
% get rhotor values
gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)';
gdat_data.rhotornorm(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.phi(end,it)));
% get rhovol values
gdat_data.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)';
gdat_data.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
gdat_data.rhovolnorm(:,it) = sqrt(abs(gdat_data.vol(:,it) ./ gdat_data.vol(end,it)));
gdat_data.area(:,it)=Area.value(it,2*Lpf1-1:-2:1)';
gdat_data.dareadpsi(:,it)=Area.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
gdat_data.Rmesh(:,it) = Ri.value(it,1:M_Rmesh);
gdat_data.Zmesh(:,it) = Zj.value(it,1:N_Zmesh);
gdat_data.psi2D(1:M_Rmesh,1:N_Zmesh,it) = PFM_tree.value(1:M_Rmesh,1:N_Zmesh,it);
gdat_data.pressure(:,it)=Pres.value(it,2*Lpf1-1:-2:1)';
gdat_data.dpressuredpsi(:,it)=Pres.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
if ~isempty(Jpol.value)
gdat_data.jpol(:,it)=Jpol.value(it,2*Lpf1-1:-2:1)';
gdat_data.djpolpsi(:,it)=Jpol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
else
gdat_data.jpol = [];
gdat_data.djpolpsi = [];
end
if ~strcmp(DIAG,'TRE')
gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)';
else
gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)'/2/pi;
end
gdat_data.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1);
gdat_data.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1);
gdat_data.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1);
if ~isempty(Rinv.value)
gdat_data.rinv(:,it) = Rinv.value(it,Lpf1:-1:1)';
else
gdat_data.rinv = [];
end
if ~isempty(R2inv.value)
gdat_data.r2inv(:,it) = R2inv.value(it,Lpf1:-1:1)';
else
gdat_data.r2inv = [];
end
if ~isempty(Bave.value)
gdat_data.bave(:,it) = Bave.value(it,Lpf1:-1:1)';
else
gdat_data.bave = [];
end
if ~isempty(B2ave.value)
gdat_data.b2ave(:,it) = B2ave.value(it,Lpf1:-1:1)';
else
gdat_data.b2ave = [];
end
if ~isempty(FTRA.value)
gdat_data.ftra(:,it) = FTRA.value(it,Lpf1:-1:1)';
else
gdat_data.ftra = [];
end
%
end
gdat_data.x = gdat_data.rhopolnorm;
%
try
[equil_Rcoil,e]=rdaAUG_eff(shot,DIAG,'Rcl',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.Rcoils = equil_Rcoil.value;
[equil_Zcoil,e]=rdaAUG_eff(shot,DIAG,'Zcl',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.Zcoils = equil_Zcoil.value;
catch
gdat_data.Rcoils = [];
gdat_data.Zcoils = [];
end
%
if strcmp(DIAG,'IDE')
[IpiPSI,e]=rdaAUG_eff(shot,'IDG','Itor',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
if length(IpiPSI.value)~=NTIME
disp('Problem with nb time points of IDG/Itor with respect to IDE NTIME, check with O. Sauter')
return
end
else
[IpiPSI,e]=rdaAUG_eff(shot,DIAG,'IpiPSI',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
end
gdat_data.Ip = IpiPSI.value(1:NTIME);
%
gdat_data.data = gdat_data.qvalue; % put q in data
gdat_data.units=[]; % not applicable
gdat_data.data_fullpath = [DIAG ' from exp_name: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
gdat_data.cocos = 17; % should check FPP
gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits{1} = 'rhopolnorm';
gdat_data.dimunits{2} = 'time [s]';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ids'}
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 = [];
disp('need an ids name in ''source'' parameter');
disp('check substructure gdat_params.source_available for an ids list');
end
ids_gen_ok = exist('ids_gen');
equil_empty = struct([]);
if ids_gen_ok ~= 2
ids_struct_saved = 'ids_structures_20190312.mat';
if ~exist(ids_struct_saved,'file')
warning(['function ids_gen not available neither file ids_structures_20190312.mat thus cannot create empty ids: ids_gen(''' ids_top_name ''')']);
return
else
eval(['load ' ids_struct_saved ])
if isfield(ids_structures,ids_top_name)
equil_empty = ids_structures.(ids_top_name);
else
if ~isempty(ids_top_name);
warning(['ids ''' ids_top_name ''' does not exist in ids_structures saved in ' ids_struct_saved]);
end
gdat_data.gdat_params.source_available = ids_list;
return
end
end
else
equil_empty = ids_gen(ids_top_name);
end
try
if ~isempty(shot)
eval(['[ids_top,ids_top_description]=aug_get_ids_' ids_top_name '(shot,equil_empty);'])
gdat_data.(ids_top_name) = ids_top;
gdat_data.([ids_top_name '_description']) = ids_top_description;
else
gdat_data.(ids_top_name) = equil_empty;
gdat_data.([ids_top_name '_description']) = ['shot empty so return default empty structure for ' ids_top_name];
end
catch ME_aug_get_ids
getReport(ME_aug_get_ids)
disp(['there is a problem with: aug_get_ids_' ids_top_name ...
' , may be check if it exists in your path or test it by itself'])
gdat_data.(ids_top_name) = equil_empty;
gdat_data.([ids_top_name '_description']) = getReport(ME_aug_get_ids);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ne','te'}
exp_name_eff = gdat_data.gdat_params.exp_name;
% ne or Te from Thomson data on raw z mesh vs (z,t)
nodenameeff = [upper(data_request_eff(1)) 'e_c'];
node_child_nameeff = [upper(data_request_eff(1)) 'e_core'];
[a,error_status]=rdaAUG_eff(shot,'VTA',nodenameeff,exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
if isempty(a.data) || isempty(a.t) || error_status>0
if gdat_params.nverbose>=3;
a
disp(['with data_request= ' data_request_eff])
end
return
end
gdat_data.(lower(node_child_nameeff)).data = a.data;
gdat_data.(lower(node_child_nameeff)).t = a.t;
gdat_data.t = a.t;
if any(strcmp(fieldnames(a),'units'))
gdat_data.units=a.units;
end
[alow,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'elow_c'],exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
[aup,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'eupp_c'],exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.(lower(node_child_nameeff)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff)).data,gdat_data.(lower(node_child_nameeff)).data-alow.value);
[a,error_status]=rdaAUG_eff(shot,'VTA','R_core',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.(lower(node_child_nameeff)).r = repmat(a.data,size(gdat_data.(lower(node_child_nameeff)).data,1),1);
[a,error_status]=rdaAUG_eff(shot,'VTA','Z_core',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.(lower(node_child_nameeff)).z = repmat(a.data',1,size(gdat_data.(lower(node_child_nameeff)).data,2));
gdat_data.dimunits=[{'Z [m]'} ; {'time [s]'}];
gdat_data.data_fullpath=[data_request_eff ' from VTA/' upper(data_request_eff(1)) 'e_c and ' upper(data_request_eff(1)) 'e_e'];
nb_core = size(gdat_data.(lower(node_child_nameeff)).z,1);
gdat_data.data = [];
gdat_data.data(1:nb_core,:) = gdat_data.(lower(node_child_nameeff)).data(1:nb_core,:);
gdat_data.dim=[{gdat_data.(lower(node_child_nameeff)).z};{gdat_data.(lower(node_child_nameeff)).t}];
gdat_data.error_bar(1:nb_core,:) = gdat_data.(lower(node_child_nameeff)).error_bar(1:nb_core,:);
% add edge part
nodenameeff_e = [upper(data_request_eff(1)) 'e_e'];
node_child_nameeff_e = [upper(data_request_eff(1)) 'e_edge'];
[a,error_status]=rdaAUG_eff(shot,'VTA',nodenameeff_e,exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.(lower(node_child_nameeff_e)).data = a.data;
ij_edge_notok = find(a.data>5e21);
gdat_data.(lower(node_child_nameeff_e)).data(ij_edge_notok) = NaN;
gdat_data.(lower(node_child_nameeff_e)).t = a.t;
if ~isempty(a.data)
[a,error_status]=rdaAUG_eff(shot,'VTA','R_edge',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.(lower(node_child_nameeff_e)).r = repmat(a.data,size(gdat_data.(lower(node_child_nameeff_e)).data,1),1);
[a,error_status]=rdaAUG_eff(shot,'VTA','Z_edge',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.(lower(node_child_nameeff_e)).z = repmat(a.data',1,size(gdat_data.(lower(node_child_nameeff_e)).data,2));
nb_edge = size(gdat_data.(lower(node_child_nameeff_e)).z,1);
iaaa=iround_os(gdat_data.(lower(node_child_nameeff_e)).t,gdat_data.(lower(node_child_nameeff)).t);
gdat_data.data(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).data(1:nb_edge,iaaa);
gdat_data.dim{1}(nb_core+1:nb_core+nb_edge,:)=gdat_data.(lower(node_child_nameeff_e)).z(1:nb_edge,iaaa);
[alow,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'elow_e'],exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
[aup,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'eupp_e'],exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.(lower(node_child_nameeff_e)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff_e)).data, ...
gdat_data.(lower(node_child_nameeff_e)).data-alow.value);
gdat_data.error_bar(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).error_bar(1:nb_edge,iaaa);
else
gdat_data.(lower(node_child_nameeff_e)).data = [];
gdat_data.(lower(node_child_nameeff_e)).t = [];
gdat_data.(lower(node_child_nameeff_e)).error_bar = [];
gdat_data.(lower(node_child_nameeff_e)).r = [];
gdat_data.(lower(node_child_nameeff_e)).z = [];
end
gdat_data.x=gdat_data.dim{1};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ne_rho', 'te_rho', 'nete_rho'}
sources_available = {'VTA','IDA','TRA'}; % TRA to be included
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = 'VTA'; % default source
end
sources_exp_name_available = 'AUGD';
if isfield(gdat_data.gdat_params,'source_exp_name') && ~isempty(gdat_data.gdat_params.source_exp_name)
sources_exp_name_available = gdat_data.gdat_params.source_exp_name;
end
gdat_data.gdat_params.source_exp_name = sources_exp_name_available;
gdat_data.gdat_params.source = upper(gdat_data.gdat_params.source);
if strcmp(gdat_data.gdat_params.source(1:2),'EQ');
warning('set equilibrium choice in ''equil'' parameter and not source, moved to equil')
gdat_data.gdat_params.equil = gdat_data.gdat_params.source;
gdat_data.gdat_params.source = [];
end
if ~isempty(intersect({'THOMSON'},gdat_data.gdat_params.source)); gdat_data.gdat_params.source = 'VTA'; end
if ~isempty(intersect({'TRANSP'},gdat_data.gdat_params.source)); gdat_data.gdat_params.source = 'TRA'; end
if ~isempty(intersect({'IDE'},gdat_data.gdat_params.source)); gdat_data.gdat_params.source = 'IDA'; end
if isempty(intersect(sources_available,gdat_data.gdat_params.source))
error(['available source choices: ' sources_available{:} char(10) ' source: ' gdat_data.gdat_params.source ...
' not implemented, ask O. Sauter' char(10)]);
elseif ~strcmp(gdat_data.gdat_params.source,'VTA')
disp(['At this stage .te, .ne, .fit refer to VTA projected on given equil choice. ' char(10) ...
'IDA or TRA profiles are added in .ida or .transp respectively']);
end
if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit) || ~isnumeric(gdat_data.gdat_params.fit)
if strcmp(gdat_data.gdat_params.source,'VTA')
gdat_data.gdat_params.fit = 1; % default do fit (only with VTA source requested)
else
% fit profiles filled in by IDA, TRA, ... profiles
gdat_data.gdat_params.fit = 0;
end
end
%
params_eff = gdat_data.gdat_params;
params_eff.data_request=data_request_eff(1:2); % start with ne if nete_rho
% get raw data
[gdat_data,params_kin,error_status]=gdat_aug(shot,params_eff);
gdat_data.gdat_params.data_request=data_request_eff;
gdat_data.gdat_request=data_request_eff;
if error_status>0
if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
return
end
% add rho coordinates
params_eff.data_request='equil';
[equil,params_equil,error_status]=gdat_aug(shot,params_eff);
if error_status>0
if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
return
end
gdat_data.gdat_params.equil = params_equil.equil;
%gdat_data.equil = equil;
% core
node_child_nameeff = [upper(data_request_eff(1)) 'e_core'];
inb_chord_core=size(gdat_data.(lower(node_child_nameeff)).r,1);
inb_time_core=size(gdat_data.(lower(node_child_nameeff)).r,2);
psi_out_core = NaN*ones(inb_chord_core,inb_time_core);
rhopolnorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
rhotornorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
rhovolnorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
% edge
node_child_nameeff_e = [upper(data_request_eff(1)) 'e_edge'];
inb_chord_edge=size(gdat_data.(lower(node_child_nameeff_e)).r,1);
inb_time_edge=size(gdat_data.(lower(node_child_nameeff_e)).r,2);
psi_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
rhopolnorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
rhotornorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
rhovolnorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[min(gdat_data.(lower(node_child_nameeff)).t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.(lower(node_child_nameeff)).t(end)+0.1)];
for itequil=1:length(time_equil)-1
rr=equil.Rmesh(:,itequil);
zz=equil.Zmesh(:,itequil);
psirz_in = equil.psi2D(:,:,itequil);
it_core_inequil = find(gdat_data.(lower(node_child_nameeff)).t>time_equil(itequil) & gdat_data.(lower(node_child_nameeff)).t<=time_equil(itequil+1));
if ~isempty(it_core_inequil)
rout_core=gdat_data.(lower(node_child_nameeff)).r(:,it_core_inequil);
zout_core=gdat_data.(lower(node_child_nameeff)).z(:,it_core_inequil);
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core);
psi_out_core(:,it_core_inequil) = reshape(psi_at_routzout,inb_chord_core,length(it_core_inequil));
rhopolnorm_out_core(:,it_core_inequil) = sqrt((psi_out_core(:,it_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
for it_cx=1:length(it_core_inequil)
rhotornorm_out_core(:,it_core_inequil(it_cx)) = ...
interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
rhovolnorm_out_core(:,it_core_inequil(it_cx)) = ...
interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
% edge
it_edge_inequil = find(gdat_data.(lower(node_child_nameeff_e)).t>time_equil(itequil) & gdat_data.(lower(node_child_nameeff_e)).t<=time_equil(itequil+1));
if ~isempty(it_edge_inequil)
rout_edge=gdat_data.(lower(node_child_nameeff_e)).r(:,it_edge_inequil);
zout_edge=gdat_data.(lower(node_child_nameeff_e)).z(:,it_edge_inequil);
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_edge,zout_edge);
psi_out_edge(:,it_edge_inequil) = reshape(psi_at_routzout,inb_chord_edge,length(it_edge_inequil));
rhopolnorm_out_edge(:,it_edge_inequil) = sqrt((psi_out_edge(:,it_edge_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
for it_cx=1:length(it_edge_inequil)
rhotornorm_out_edge(:,it_edge_inequil(it_cx)) = ...
interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
rhovolnorm_out_edge(:,it_edge_inequil(it_cx)) = ...
interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
end
gdat_data.(lower(node_child_nameeff)).psi = psi_out_core;
gdat_data.(lower(node_child_nameeff)).rhopolnorm = rhopolnorm_out_core;
gdat_data.(lower(node_child_nameeff)).rhotornorm = rhotornorm_out_core;
gdat_data.(lower(node_child_nameeff)).rhovolnorm = rhovolnorm_out_core;
gdat_data.(lower(node_child_nameeff_e)).psi = psi_out_edge;
gdat_data.(lower(node_child_nameeff_e)).rhopolnorm = rhopolnorm_out_edge;
gdat_data.(lower(node_child_nameeff_e)).rhotornorm = rhotornorm_out_edge;
gdat_data.(lower(node_child_nameeff_e)).rhovolnorm = rhovolnorm_out_edge;
% put values of rhopolnorm for dim{1} by default, all radial mesh for combined core, edge in grids_1d
gdat_data.x = gdat_data.(lower(node_child_nameeff)).rhopolnorm;
iaaa=iround_os(gdat_data.(lower(node_child_nameeff_e)).t,gdat_data.(lower(node_child_nameeff)).t);
gdat_data.x(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhopolnorm(:,iaaa);
gdat_data.dim{1} = gdat_data.x;
gdat_data.dimunits{1} = 'rhopolnorm';
gdat_data2 = get_grids_1d(gdat_data,2,1,gdat_params.nverbose);
gdat_data.grids_1d = gdat_data2.grids_1d;
% $$$ gdat_data.grids_1d.rhopolnorm = gdat_data.x;
% $$$ gdat_data.grids_1d.psi = gdat_data.(lower(node_child_nameeff)).psi;
% $$$ gdat_data.grids_1d.psi(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).psi(:,iaaa);
% $$$ gdat_data.grids_1d.rhotornorm = gdat_data.(lower(node_child_nameeff)).rhotornorm;
% $$$ gdat_data.grids_1d.rhotornorm(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhotornorm(:,iaaa);
% $$$ gdat_data.grids_1d.rhovolnorm = gdat_data.(lower(node_child_nameeff)).rhovolnorm;
% $$$ gdat_data.grids_1d.rhovolnorm(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhovolnorm(:,iaaa);
gdat_data.data_fullpath = [gdat_data.data_fullpath ' projected on equilibrium ' gdat_data.gdat_params.equil];
% if nete_rho, copy data as .ne, get .te and put pe=e ne Te in data:
gdat_data.(data_request_eff(1:2)).data = gdat_data.data;
gdat_data.(data_request_eff(1:2)).error_bar = gdat_data.error_bar;
gdat_data.(data_request_eff(1:2)).units = gdat_data.units;
gdat_data.(data_request_eff(1:2)).core = gdat_data.([data_request_eff(1:2) '_core']);
gdat_data.(data_request_eff(1:2)).edge = gdat_data.([data_request_eff(1:2) '_edge']);
gdat_data = rmfield(gdat_data,{[data_request_eff(1:2) '_core'],[data_request_eff(1:2) '_edge']});
if strcmp(data_request_eff(1:4),'nete')
params_eff.data_request=data_request_eff(3:4);
[gdat_data_te,params_kin,error_status]=gdat_aug(shot,params_eff);
gdat_data.te.data = gdat_data_te.data;
gdat_data.te.error_bar = gdat_data_te.error_bar;
gdat_data.te.units = gdat_data_te.units;
gdat_data.te.core = gdat_data_te.te_core;
gdat_data.te.edge = gdat_data_te.te_edge;
gdat_data.te.error_bar = gdat_data_te.error_bar;
gdat_data.te.core.psi = gdat_data.ne.core.psi;
gdat_data.te.core.rhopolnorm = gdat_data.ne.core.rhopolnorm;
gdat_data.te.core.rhotornorm = gdat_data.ne.core.rhotornorm;
gdat_data.te.core.rhovolnorm = gdat_data.ne.core.rhovolnorm;
gdat_data.te.edge.psi = gdat_data.ne.edge.psi;
gdat_data.te.edge.rhopolnorm = gdat_data.ne.edge.rhopolnorm;
gdat_data.te.edge.rhotornorm = gdat_data.ne.edge.rhotornorm;
gdat_data.te.edge.rhovolnorm = gdat_data.ne.edge.rhovolnorm;
% put pe in main gdat_data
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);
gdat_data.units='N/m^2; 1.6022e-19 ne Te';
gdat_data.data_fullpath=['pe=1.6e-19*ne*Te in data, .ne, .te ' gdat_data.data_fullpath(3:end)];
gdat_data.label = 'pe';
end
% defaults for fits, so user always gets std structure
gdat_data.fit.rhotornorm = []; % same for both ne and te
gdat_data.fit.rhopolnorm = [];
gdat_data.fit.t = [];
gdat_data.fit.te.data = [];
gdat_data.fit.te.drhotornorm = [];
gdat_data.fit.ne.data = [];
gdat_data.fit.ne.drhotornorm = [];
gdat_data.fit.raw.te.data = [];
gdat_data.fit.raw.te.rhotornorm = [];
gdat_data.fit.raw.ne.data = [];
gdat_data.fit.raw.ne.rhotornorm = [];
fit_tension_default = -0.1;
if isfield(gdat_data.gdat_params,'fit_tension')
fit_tension = gdat_data.gdat_params.fit_tension;
else
fit_tension = fit_tension_default;
end
if ~isstruct(fit_tension)
fit_tension_eff.te = fit_tension;
fit_tension_eff.ne = fit_tension;
fit_tension = fit_tension_eff;
else
if ~isfield(fit_tension,'te'); fit_tension.te = fit_tension_default; end
if ~isfield(fit_tension,'ne '); fit_tension.ne = fit_tension_default; end
end
gdat_data.gdat_params.fit_tension = fit_tension;
if isfield(gdat_data.gdat_params,'fit_nb_rho_points')
fit_nb_rho_points = gdat_data.gdat_params.fit_nb_rho_points;
else
fit_nb_rho_points = 201;
end
gdat_data.gdat_params.fit_nb_rho_points = fit_nb_rho_points;
%
if gdat_data.gdat_params.fit==1
% add "manual" fits
gdat_data.fit.t = gdat_data.t;
rhotornormfit = linspace(0,1,fit_nb_rho_points)';
gdat_data.fit.rhotornorm = rhotornormfit;
gdat_data.fit.rhopolnorm = NaN*ones(length(rhotornormfit),length(gdat_data.fit.t));
if any(strfind(data_request_eff,'te'))
gdat_data.fit.raw.te.data = NaN*ones(size(gdat_data.te.data));
gdat_data.fit.raw.te.error_bar = NaN*ones(size(gdat_data.te.data));
gdat_data.fit.raw.te.rhotornorm = NaN*ones(size(gdat_data.te.data));
gdat_data.fit.te.data = gdat_data.fit.rhopolnorm;
gdat_data.fit.te.drhotornorm = gdat_data.fit.rhopolnorm;
end
if any(strfind(data_request_eff,'ne'))
gdat_data.fit.raw.ne.data = NaN*ones(size(gdat_data.ne.data));
gdat_data.fit.raw.ne.error_bar = NaN*ones(size(gdat_data.ne.data));
gdat_data.fit.raw.ne.rhotornorm = NaN*ones(size(gdat_data.ne.data));
gdat_data.fit.ne.data =gdat_data.fit.rhopolnorm;
gdat_data.fit.ne.drhotornorm = gdat_data.fit.rhopolnorm;
end
for it=1:length(gdat_data.t)
% make rhotor->rhopol transformation for each time since equilibrium might have changed
gdat_data.fit.rhopolnorm(:,it)=interpos(gdat_data.grids_1d.rhotornorm_equil(:,it), ...
gdat_data.grids_1d.rhopolnorm_equil(:,it),rhotornormfit);
gdat_data.fit.rhopolnorm(:,it)=interpos(gdat_data.grids_1d.rhotornorm_equil(:,it), ...
gdat_data.grids_1d.rhopolnorm_equil(:,it),rhotornormfit);
if any(strfind(data_request_eff,'te'))
idatate = find(gdat_data.te.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05);
if length(idatate)>0
gdat_data.fit.raw.te.rhotornorm(idatate,it) = gdat_data.grids_1d.rhotornorm(idatate,it);
gdat_data.fit.raw.te.data(idatate,it) = gdat_data.te.data(idatate,it);
gdat_data.fit.raw.te.error_bar(idatate,it) = gdat_data.te.error_bar(idatate,it);
[rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhotornorm(idatate,it));
rhoeff = [0; rhoeff];
teeff = gdat_data.te.data(idatate(irhoeff),it);
te_err_eff = gdat_data.te.error_bar(idatate(irhoeff),it);
% they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
ij=find(teeff./te_err_eff>10.);
if ~isempty(ij); te_err_eff(ij) = max(max(te_err_eff),teeff(ij)./0.1); end
%
teeff = [teeff(1); teeff];
te_err_eff = [1e4; te_err_eff];
% add some points to reach 0 at 1.05
if max(rhoeff)<1.03
rhoeff = [rhoeff; 1.05];
teeff = [teeff; 0.];
te_err_eff(end) = 10.*max(te_err_eff(2:end));
te_err_eff = [te_err_eff; min(te_err_eff)]; % to give more weight for new last point
else
teeff(end) = 0.;
te_err_eff(end-1) = 10.*max(te_err_eff(2:end));
te_err_eff(end) = min(te_err_eff); % to give more weight for new value of last point
end
[gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhotornorm(:,it)] = interpos(rhoeff,teeff,rhotornormfit,fit_tension.te,[1 0],[0 0],te_err_eff);
% avoid neg points or positive edge gradient
if ~isempty(find(gdat_data.fit.te.data(:,it)<0)) || gdat_data.fit.te.drhotornorm(end,it) > 0
ij= find(rhoeff>=0.85 & rhoeff<1.04);
if ~isempty(ij)
te_err_eff(ij) = 100.*min(te_err_eff);
else
te_err_eff(end-min(5,length(te_err_eff)-3):end-1) = 100.*min(te_err_eff);
end
[gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhotornorm(:,it)] = interpos(rhoeff,teeff,rhotornormfit,fit_tension.te,[1 0],[0 0],te_err_eff);
end
end
end
if any(strfind(data_request_eff,'ne'))
idatane = find(gdat_data.ne.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05);
if length(idatane)>0
gdat_data.fit.raw.ne.rhotornorm(idatane,it) = gdat_data.grids_1d.rhotornorm(idatane,it);
gdat_data.fit.raw.ne.data(idatane,it) = gdat_data.ne.data(idatane,it);
gdat_data.fit.raw.ne.error_bar(idatane,it) = gdat_data.ne.error_bar(idatane,it);
[rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhotornorm(idatane,it));
rhoeff = [0; rhoeff];
% they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
neeff = gdat_data.ne.data(idatane(irhoeff),it);
ne_err_eff = gdat_data.ne.error_bar(idatane(irhoeff),it);
ij=find(neeff./ne_err_eff>10.);
if ~isempty(ij); ne_err_eff(ij) = max(max(ne_err_eff),neeff(ij)./0.1); end
%
neeff = [neeff(1); neeff];
ne_err_eff = [1e21; ne_err_eff];
% add some points to reach 0 at 1.05
if max(rhoeff)<1.03
rhoeff = [rhoeff; 1.05];
neeff = [neeff; 0.];
ne_err_eff(end) = 10.*max(ne_err_eff(2:end)); % to give more weight for new last point
ne_err_eff = [ne_err_eff; min(ne_err_eff)];
else
neeff(end) = 0.;
ne_err_eff(end-1) = 10.*max(ne_err_eff(2:end));
ne_err_eff(end) = min(ne_err_eff); % to give more weight for new value of last point
end
[gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhotornorm(:,it)] = interpos(rhoeff,neeff,rhotornormfit,fit_tension.ne,[1 0],[0 0],ne_err_eff);
% avoid neg points or positive edge gradient
if ~isempty(find(gdat_data.fit.ne.data(:,it)<0)) || gdat_data.fit.ne.drhotornorm(end,it) > 0
ij= find(rhoeff>=0.85 & rhoeff<1.04);
if ~isempty(ij)
ne_err_eff(ij) = 100.*min(ne_err_eff);
else
ne_err_eff(end-min(5,length(ne_err_eff)-3):end-1) = 100.*min(ne_err_eff);
end
[gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhotornorm(:,it)] = interpos(rhoeff,neeff,rhotornormfit,fit_tension.ne,[1 0],[0 0],ne_err_eff);
end
end
end
end
end
% add other profiles in fit (without raw data part)
switch gdat_data.gdat_params.source
case 'IDA'
gdat_data.fit.full_path = 'from IDA: Te, Te_unc, ne, ne_unc, rhop and rhot';
params_eff.data_request={'IDA','','','area-base:rhop'};
[gdat_data_ida_rhop,params_kin,error_status]=gdat_aug(shot,params_eff);
gdat_data.fit.rhopolnorm = gdat_data_ida_rhop.data;
params_eff.data_request={'IDA','rhot'}; % it is actually rhotor_norm
[gdat_data_te_ida_rhot,params_kin,error_status]=gdat_aug(shot,params_eff);
gdat_data.fit.rhotornorm = gdat_data_te_ida_rhot.data;
if strcmp(data_request_eff(1:4),'nete') || strcmp(data_request_eff(1:2),'te')
% Te
params_eff.data_request={'IDA','Te'};
[gdat_data_te_ida,params_kin,error_status]=gdat_aug(shot,params_eff);
gdat_data.fit.te.data = gdat_data_te_ida.data;
gdat_data.fit.t = gdat_data_te_ida.t;
params_eff.data_request={'IDA','Te_unc'};
[gdat_data_te_ida_err,params_kin,error_status]=gdat_aug(shot,params_eff);
gdat_data.fit.te.error_bar = gdat_data_te_ida_err.data;
end
if strcmp(data_request_eff(1:4),'nete') || strcmp(data_request_eff(1:2),'ne')
% ne
params_eff.data_request={'IDA','ne'};
[gdat_data_ne_ida,params_kin,error_status]=gdat_aug(shot,params_eff);
gdat_data.fit.ne.data = gdat_data_ne_ida.data;
gdat_data.fit.t = gdat_data_ne_ida.t;
params_eff.data_request={'IDA','ne_unc'};
[gdat_data_ne_ida_err,params_kin,error_status]=gdat_aug(shot,params_eff);
gdat_data.fit.ne.error_bar = gdat_data_ne_ida_err.data; % time as 2nd dim
end
case 'TRA'
gdat_data.fit.full_path = 'from TRA: TE, NE, rhop from PLFLX and XB (.x is rhotornorm mid-points, XB end-point)';
% (tra.XB.data(1:end,it)'+[0 ,tra.XB.data(1:end-1,it)'])/2 = tra.TE.x(1:end,it)'
if strcmp(data_request_eff(1:4),'nete') || strcmp(data_request_eff(1:2),'te')
% Te
params_eff.data_request={'TRA','TE',gdat_data.gdat_params.source_exp_name};
[gdat_data_te_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
gdat_data.fit.te.data = gdat_data_te_tra.data;
gdat_data.fit.t = gdat_data_te_tra.t;
gdat_data.fit.rhotornorm = gdat_data_te_tra.x;
params_eff.data_request={'TRA','PLFLX',gdat_data.gdat_params.source_exp_name};
[gdat_data_plflx_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
rhopol_endpoints = sqrt((gdat_data_plflx_tra.data-0)./(repmat(gdat_data_plflx_tra.data(end,:),size(gdat_data_plflx_tra.data,1),1)-0.));
for it=1:size(gdat_data_plflx_tra.data,2)
gdat_data.fit.rhopolnorm(:,it) = interpos([0. ; gdat_data_plflx_tra.x(:,it)],[0. ; rhopol_endpoints(:,it)],gdat_data.fit.rhotornorm(:,it));
end
end
if strcmp(data_request_eff(1:4),'nete') || strcmp(data_request_eff(1:2),'ne')
% ne
params_eff.data_request={'TRA','NE',gdat_data.gdat_params.source_exp_name};
[gdat_data_ne_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
gdat_data.fit.ne.data = gdat_data_ne_tra.data .* 1e6; % ne in cm^3 !!
if ~strcmp(data_request_eff(1:4),'nete')
gdat_data.fit.t = gdat_data.t;
gdat_data.fit.rhotornorm = gdat_data_ne_tra.x;
params_eff.data_request={'TRA','PLFLX',gdat_data.gdat_params.source_exp_name};
[gdat_data_plflx_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
rhopol_endpoints = sqrt((gdat_data_plflx_tra.data-0)./(repmat(gdat_data_plflx_tra.data(end,:),size(gdat_data_plflx_tra.data,1),1)-0.));
for it=1:size(gdat_data_plflx_tra.data,2)
gdat_data.fit.rhopolnorm(:,it) = interpos([0. ; gdat_data_plflx_tra.x(:,it)],[0. ; rhopol_endpoints(:,it)],gdat_data.fit.rhotornorm(:,it));
end
end
end
otherwise
if ~strcmp(gdat_data.gdat_params.source,'VTA')
disp(['profiles ''fit'' from ' gdat_data.gdat_params.source ' not defined. At this stage only: ' sources_available{:} ...
' are available. Ask O. Sauter if need be']);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'pgyro'}
extra_arg_sf2sig_eff_string = '';
if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
extra_arg_sf2sig_eff_string = [',' gdat_data.gdat_params.extra_arg_sf2sig];
end
% LOAD MULTI CHANNEL DATA ECS
% powers, frequencies, etc
params_eff = gdat_data.gdat_params;
params_eff.data_request={'ECS','PECRH'};
% pgyro tot in index=9
try
gdat_data=gdat_aug(shot,params_eff);
gdat_data.data_request = data_request_eff;
gdat_data.gdat_params.data_request = data_request_eff;
catch
if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
gdat_data.data_request = data_request_eff;
return
end
nb_timepoints = length(gdat_data.t);
pgyro = NaN*ones(nb_timepoints,9);
pgyro(:,9) = reshape(gdat_data.data,nb_timepoints,1);
gdat_data.data = pgyro;
labels{9} = 'ECtot';
for i=1:4
% "old" ECRH1 gyrotrons: gyro 1 to 4 in pgyro
params_eff.data_request={'ECS',['PG' num2str(i)]};
gdat_data_i=gdat_aug(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i} = gdat_data_i;
gdat_data.data(:,i) = reshape(gdat_data_i.data,nb_timepoints,1);
gdat_data.dim = [{gdat_data_i.t} {[1:9]}];
if max(gdat_data_i.data) > 0.
labels{i} = ['EC_' num2str(i)];
end
end
sys1_source = 'P_sy1_g';
if shot > 33725; sys1_source = 'P_sy3_g'; end
try
[a,e]=rdaAUG_eff(shot,'ECS',[sys1_source num2str(i)],gdat_data.gdat_params.exp_name,[], ...
gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'gyr_freq']);
% eval(['a = sf2par(''ECS'',shot,''gyr_freq'','sys1_source num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
catch
% gyr_freq not present (old shots for example)
a=[];
end
if isempty(a)
else
gdat_data.ec{i}.freq = a;
gdat_data.freq_ec(i) = a.value;
end
try
[a,e]=rdaAUG_eff(shot,'ECS',[sys1_source num2str(i)],gdat_data.gdat_params.exp_name,[], ...
gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GPolPos']);
% eval(['a = sf2par(''ECS'',shot,''GPolPos'','sys1_source num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
catch
% GPolPos not present
a=[];
end
if isempty(a)
else
gdat_data.ec{i}.polpos = a.value;
gdat_data.polpos_ec(i) = a.value;
end
try
[a,e]=rdaAUG_eff(shot,'ECS',[sys1_source num2str(i)],gdat_data.gdat_params.exp_name,[], ...
gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GTorPos']);
% eval(['a = sf2par(''ECS'',shot,''GTorPos'','sys1_source num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
catch
a=[];
end
if isempty(a)
else
gdat_data.ec{i}.torpos = a.value;
gdat_data.torpos_ec(i) = a.value;
end
% "new" ECRH2 gyrotrons: gyro 5 to 8 in pgyro
params_eff.data_request={'ECS',['PG' num2str(i) 'N']};
gdat_data_i=gdat_aug(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i+4} = gdat_data_i;
gdat_data.data(:,i+4) = reshape(gdat_data_i.data,nb_timepoints,1);
gdat_data.dim = [{gdat_data_i.t} {[1:9]}];
if max(gdat_data_i.data) > 0.
labels{i+4} = ['EC_' num2str(i+4)];
end
end
try
[a,e]=rdaAUG_eff(shot,'ECS',['P_sy2_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ...
gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'gyr_freq']);
% eval(['a = sf2par(''ECS'',shot,''gyr_freq'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
catch
a=[];
end
if isempty(a)
else
gdat_data.ec{i+4}.freq = a;
gdat_data.freq_ec(i+4) = a.value;
end
try
[a,e]=rdaAUG_eff(shot,'ECS',['P_sy2_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ...
gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GPolPos']);
% eval(['a = sf2par(''ECS'',shot,''GPolPos'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
catch
a=[];
end
if isempty(a)
else
gdat_data.ec{i+4}.polpos = a.value;
gdat_data.polpos_ec(i+4) = a.value;
end
try
[a,e]=rdaAUG_eff(shot,'ECS',['P_sy2_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ...
gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GTorPos']);
% eval(['a = sf2par(''ECS'',shot,''GTorPos'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']);
catch
a=[];
end
if isempty(a)
else
gdat_data.ec{i+4}.torpos = a.value;
gdat_data.torpos_ec(i+4) = a.value;
end
params_eff.data_request={'ECN',['G' num2str(i) 'POL']};
gdat_data_i=gdat_aug(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i+4}.gpol_ec = gdat_data_i;
end
params_eff.data_request={'ECN',['G' num2str(i) 'TOR']};
gdat_data_i=gdat_aug(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i+4}.gtor_ec = gdat_data_i;
end
params_eff.data_request={'ECN',['G' num2str(i) 'PO4']};
gdat_data_i=gdat_aug(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i+4}.gpo4_ec = gdat_data_i;
end
params_eff.data_request={'ECN',['G' num2str(i) 'PO8']};
gdat_data_i=gdat_aug(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i+4}.gpo8_ec = gdat_data_i;
end
end
if ~isempty(gdat_data.dim)
gdat_data.t = gdat_data.dim{1};
gdat_data.x = gdat_data.dim{2};
gdat_data.dimunits=[{'time [s]'} {'ECRH1(1:4) ECRH2(1:4) ECtot'}];
gdat_data.units='W';
gdat_data.freq_ech_units = 'GHz';
gdat_data.data_fullpath=['ECS/' 'PGi and PGiN, etc'];
icount=0;
for i=1:length(labels)
if ~isempty(labels{i})
icount=icount+1;
gdat_data.label{icount} = labels{i};
end
end
else
gdat_data.freq_ech_units =[]';
gdat_data.ec = [];
gdat_data.freq_ec = [];
gdat_data.polpos_ec = [];
gdat_data.torpos_ec = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'powers'}
sources_avail = {'ohm','ec','nbi','ic','rad'}; % note should allow ech, nbh, ohmic in parameter sources
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
fields_to_copy = {'data','units','dim','dimunits','t','x','data_fullpath','label','help','gdat_params'};
fields_to_not_copy = {'shot','gdat_request'};
% total of each source in .data, but full data in subfield like pgyro in .ec, to check for nbi
params_eff = gdat_data.gdat_params;
% ohmic, use its time-base
params_eff.data_request={'TOT','P_OH'};
try
ohm=gdat_aug(shot,params_eff);
catch
ohm.data = [];
ohm.dim = [];
end
if ~isempty(ohm.data) && ~isempty(ohm.dim)
for i=1:length(fields_to_copy)
if isfield(ohm,fields_to_copy{i})
gdat_data.ohm.(fields_to_copy{i}) = ohm.(fields_to_copy{i});
end
end
gdat_data.ohm.raw_data = gdat_data.ohm.data;
else
if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
return
end
mapping_for_aug.timedim = 1; mapping_for_aug.gdat_timedim = 1;
taus = -10;
%
% add each source in main.data, on ohm time array
gdat_data.t = gdat_data.ohm.t;
gdat_data.dim{1} = gdat_data.t;
gdat_data.dimunits{1} = 's';
gdat_data.ohm.data = interpos(gdat_data.t,gdat_data.ohm.raw_data,5.*taus);
gdat_data.data = reshape(gdat_data.ohm.data,length(gdat_data.t),1);
gdat_data.ohm.tension = 5.*taus;
gdat_data.x =[1];
gdat_data.label={'P_{ohm}'};
gdat_data.units = 'W';
%
if any(strmatch('ec',gdat_data.gdat_params.source))
% ec
params_eff.data_request={'ECS','PECRH'};
params_eff.data_request='pgyro';
try
ec=gdat_aug(shot,params_eff);
catch
end
if ~isempty(ec.data) && ~isempty(ec.dim)
for i=1:length(fields_to_copy)
% if has pgyro, use not_copy
if isfield(ec,fields_to_copy{i}) && ~any(strmatch(fields_to_not_copy,fields_to_copy{i}))
gdat_data.ec.(fields_to_copy{i}) = ec.(fields_to_copy{i});
end
end
gdat_data.data(:,end+1) = interpos(-21,gdat_data.ec.t,gdat_data.ec.data(:,end),gdat_data.t);
gdat_data.x(end+1) =gdat_data.x(end)+1;
gdat_data.label{end+1}='P_{ec}';
end
end
%
if any(strmatch('nb',gdat_data.gdat_params.source))
% nbi
params_eff.data_request={'NIS','PNIQ'};
try
nbi=gdat_aug(shot,params_eff);
catch
nbi.data = [];
nbi.dim = [];
end
if ~isempty(nbi.data) && ~isempty(nbi.dim)
for i=1:length(fields_to_copy)
if isfield(nbi,fields_to_copy{i})
gdat_data.nbi.(fields_to_copy{i}) = nbi.(fields_to_copy{i});
end
end
% add to main with linear interpolation and 0 for extrapolated values
gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi.t,gdat_data.nbi.data(:,end),gdat_data.t);
gdat_data.x(end+1) =gdat_data.x(end)+1;
gdat_data.label{end+1}='P_{nbi}';
end
end
%
if any(strmatch('ic',gdat_data.gdat_params.source))
% ic
params_eff.data_request={'ICP','PICRN'};
try
ic=gdat_aug(shot,params_eff);
catch
ic.data = [];
ic.dim = [];
end
if ~isempty(ic.data) && ~isempty(ic.dim)
for i=1:length(fields_to_copy)
if isfield(ic,fields_to_copy{i})
gdat_data.ic.(fields_to_copy{i}) = ic.(fields_to_copy{i});
end
end
gdat_data.data(:,end+1) = interpos(-21,gdat_data.ic.t,gdat_data.ic.data,gdat_data.t);
gdat_data.x(end+1) =gdat_data.x(end)+1;
gdat_data.label{end+1}='P_{ic}';
end
end
index_prad = [];
if any(strmatch('rad',gdat_data.gdat_params.source))
% radiated power
params_eff.data_request='prad';
try
rad=gdat_aug(shot,params_eff);
catch
end
if ~isempty(rad.data) && ~isempty(rad.dim)
for i=1:length(fields_to_copy)
if isfield(rad,fields_to_copy{i})
gdat_data.rad.(fields_to_copy{i}) = rad.(fields_to_copy{i});
end
end
% 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,gdat_data.t);
index_prad = size(gdat_data.data,2); % to not add for ptot
gdat_data.x(end+1) =gdat_data.x(end)+1;
gdat_data.label{end+1}='P_{rad}';
end
end
% add tot power
index_noprad = setdiff([1:size(gdat_data.data,2)],index_prad);
gdat_data.data(:,end+1) = sum(gdat_data.data(:,index_noprad),2);
gdat_data.label{end+1}='P_{tot,h}';
gdat_data.x(end+1) =gdat_data.x(end)+1;
gdat_data.dim{2} = gdat_data.x; gdat_data.dimunits{2} = '';
gdat_data.data_fullpath = 'tot powers from each sources, and total heating power in .data(:,end)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'q_rho'}
[gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data);
% since Lpf depends on time, need to load all first and then loop over time for easier mapping
[qpsi,e]=rdaAUG_eff(shot,DIAG,'Qpsi',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
ndimrho = size(qpsi.data,2);
if ndimrho==NTIME_Lpf
% data seems to be transposed
ndimrho = size(qpsi.data,1);
itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
else
itotransposeback = 0;
end
qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback);
ijnan=find(isnan(qpsi.value));
qpsi.value(ijnan)=0;
[psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
[phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
[Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
% seems "LCFS" q-value is far too large, limit to some max (when diverted)
max_qValue = 30.; % Note could just put a NaN on LCFS value since ill-defined when diverted
for it=1:NTIME
Lpf1 = Lpf1_t(it);
% Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
% change it to (radial,time) and use only Lpf+1 points up to LCFS
ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially
% set NaNs to zeroes
if qpsi.value(ijok(1),1)<0
gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
else
gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
end
% get x values
psi_it=psi_tree.value(it,Lpf1:-1:1)';
gdat_data.psi_axis(it)= psi_it(1);
gdat_data.psi_lcfs(it)= psi_it(end);
gdat_data.rhopolnorm(:,it) = sqrt(abs((psi_it-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it))));
if strcmp(DIAG,'EQR');
% q value has only a few values and from center to edge, assume they are from central rhopol values on
% But they are every other point starting from 3rd
ijk=find(gdat_data.qvalue(:,it)~=0);
if length(ijk)>2
% now shots have non-zero axis values in eqr
rhoeff=gdat_data.rhopolnorm(ijk,it);
qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above
if ijk(1)>1
rhoeff = [0.; rhoeff];
qeff = [qeff(1) ;qeff];
end
ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it)));
qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]);
else
qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
end
gdat_data.qvalue(:,it) = qfit;
end
% get rhotor values
phi_it = phi_tree.value(it,Lpf1:-1:1)';
gdat_data.rhotornorm(:,it) = sqrt(abs(phi_it ./ phi_it(end)));
% get rhovol values
vol_it=Vol.value(it,2*Lpf1-1:-2:1)';
gdat_data.rhovolnorm(:,it) = sqrt(abs(vol_it ./ vol_it(end)));
% Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
% change it to (radial,time) and use only Lpf+1 points up to LCFS
ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially
% max value 1e6, change to 2*extrapolation (was max_qValue before)
q_edge=interpos(gdat_data.rhotornorm(1:end-1,it),qpsi.value(it,Lpf1:-1:2),1,-0.1);
gdat_data.qvalue(:,it) = qpsi.value(it,Lpf1:-1:1)';
if abs(gdat_data.qvalue(end,it)) > 1e3
% assume diverted
gdat_data.qvalue(end,it) = 2. * q_edge;
end
% $$$ if qpsi.value(ijok(1),1)<0
% $$$ gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
% $$$ else
% $$$ gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
% $$$ end
end
gdat_data.x = gdat_data.rhopolnorm;
% get time values
gdat_data.data = gdat_data.qvalue; % put q in data
gdat_data.units=[]; % not applicable
gdat_data.data_fullpath = [DIAG ' from exp_name: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
gdat_data.cocos = 17; % should check FPP
gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits{1} = 'rhopolnorm';
gdat_data.dimunits{2} = 'time [s]';
% add grids_1d to have rhotor, etc
gdat_data = get_grids_1d(gdat_data,2,1,gdat_params.nverbose);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'psi_axis', 'psi_edge'}
if strcmp(upper(gdat_data.gdat_params.equil),'FPG')
gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev;
gdat_params_eff.data_request = [{'FPG'},{'fax-bnd'},{'AUGD'}];
gdat_data = gdat_aug(gdat_data.shot,gdat_params_eff);
gdat_data.gdat_params = gdat_params_prev;
gdat_data.label = 'psi\_axis-psi\_edge';
gdat_data.data_fullpath = [gdat_data.data_fullpath ' yields only psi\_axis-psi\_edge from FPG/fax-bnd'];
else
[gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data);
% since Lpf depends on time, need to load all first and then loop over time for easier mapping
[psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
ndimrho = size(psi_tree.data,2);
if ndimrho==NTIME_Lpf
% data seems to be transposed
ndimrho = size(psi_tree.data,1);
itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
else
itotransposeback = 0;
end
psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
ijnan=find(isnan(psi_tree.value));
psi_tree.value(ijnan)=0;
gdat_data.dim{1} = gdat_data.t;
gdat_data.dimunits{1} = 'time [s]';
for it=1:NTIME
Lpf1 =min(Lpf1_t(it),size(psi_tree.value,2));
psi_it=psi_tree.value(it,Lpf1:-1:1)';
if strcmp(data_request_eff,'psi_axis')
gdat_data.data(it)= psi_it(1);
elseif strcmp(data_request_eff,'psi_edge')
gdat_data.data(it)= psi_it(end);
else
end
end
gdat_data.units = psi_tree.units;
gdat_data.data_fullpath = [DIAG '/PFL extract ' data_request_eff];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'raptor'}
sources_avail = {'observer','predictive'};
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
rap_signals={'ECE_f','Hmode','Ip','P_LH','Phib','beta','cs','exflag','g_f','it','li3','ne_f','res','tcomp','ECEm', ...
'F','Vp','chie','dk_k','dpsi','dte','g1','g23or','g2','g3','iota','jbs','jec','jnb','jpar','ne','ni', ...
'p','paux','pbrem','pec','pei','pnbe','poh','prad','psi','q','rhtECE','sdk_k','shear','signeo','sxk_k', ...
'te','ti','upl','xk_k','ze'};
rap_desc_signals={'rts:Dia/RAPTOR/ECEflag.val','rts:Dia/RAPTOR/Hmode.val','rts:Dia/RAPTOR/Ip.val','rts:Dia/RAPTOR/P_LH.val',...
'rts:Dia/RAPTOR/Phib.val','rts:Dia/RAPTOR/beta.val','rts:Dia/RAPTOR/conf_state.val', ...
'rts:Dia/RAPTOR/exitflag.val','rts:Dia/RAPTOR/gflag.val','rts:Dia/RAPTOR/it.val', ...
'rts:Dia/RAPTOR/li3.val','rts:Dia/RAPTOR/neflag.val','rts:Dia/RAPTOR/res.val', ...
'rts:Dia/RAPTOR/dtcomp.val','rts:Dia/RAPTOR/ECEmask.val','rts:Dia/RAPTOR/F.val', ...
'rts:Dia/RAPTOR/Vp.val','rts:Dia/RAPTOR/chie.val','rts:Dia/RAPTOR/dk_k.val', ...
'rts:Dia/RAPTOR/dpsi.val','rts:Dia/RAPTOR/dte.val','rts:Dia/RAPTOR/g1.val', ...
'rts:Dia/RAPTOR/g23or.val','rts:Dia/RAPTOR/g2.val','rts:Dia/RAPTOR/g3.val', ...
'rts:Dia/RAPTOR/iota.val','rts:Dia/RAPTOR/jbs.val','rts:Dia/RAPTOR/jec.val', ...
'rts:Dia/RAPTOR/jnb.val','rts:Dia/RAPTOR/jpar.val','rts:Dia/RAPTOR/ne.val', ...
'rts:Dia/RAPTOR/ni.val','rts:Dia/RAPTOR/p.val','rts:Dia/RAPTOR/paux.val', ...
'rts:Dia/RAPTOR/pbrem.val','rts:Dia/RAPTOR/pec.val','rts:Dia/RAPTOR/pei.val', ...
'rts:Dia/RAPTOR/pnbe.val','rts:Dia/RAPTOR/poh.val','rts:Dia/RAPTOR/prad.val', ...
'rts:Dia/RAPTOR/psi.val','rts:Dia/RAPTOR/q.val','rts:Dia/RAPTOR/rhotorECE.val', ...
'rts:Dia/RAPTOR/sdk_k.val','rts:Dia/RAPTOR/shear.val','rts:Dia/RAPTOR/signeo.val', ...
'rts:Dia/RAPTOR/sxk_k.val','rts:Dia/RAPTOR/te.val','rts:Dia/RAPTOR/ti.val', ...
'rts:Dia/RAPTOR/upl.val','rts:Dia/RAPTOR/xk_k.val','rts:Dia/RAPTOR/ze.val'};
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 strcmp(gdat_data.gdat_params.source(1),'o'); gdat_data.gdat_params.source = 'observer'; end
if strcmp(gdat_data.gdat_params.source(1),'p'); gdat_data.gdat_params.source = 'predictive'; end
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 strcmp(gdat_data.gdat_params.source{i}(1),'o'); gdat_data.gdat_params.source{i} = 'observer'; end
if strcmp(gdat_data.gdat_params.source{i}(1),'p'); gdat_data.gdat_params.source{i} = 'predictive'; end
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
for i=1:length(gdat_data.gdat_params.source)
it_def=0;
for j=1:length(rap_signals)
gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev;
gdat_params_eff.data_request = {'RAP',[rap_signals{j} '_' gdat_data.gdat_params.source{i}(1)],gdat_params_eff.exp_name};
abcd = gdat_aug(gdat_data.shot,gdat_params_eff);
gdat_data.(gdat_data.gdat_params.source{i}).(rap_signals{j}) = abcd;
gdat_data.(gdat_data.gdat_params.source{i}).(rap_signals{j}).description = rap_desc_signals{j};
if strcmp(rap_signals{j},'q') && ~isempty(abcd.data)
% promote q at top level of gdat_data.(gdat_data.gdat_params.source{i})
gdat_data.(gdat_data.gdat_params.source{i}).data = abcd.data;
gdat_data.(gdat_data.gdat_params.source{i}).t = abcd.t;
gdat_data.(gdat_data.gdat_params.source{i}).x = abcd.x;
gdat_data.(gdat_data.gdat_params.source{i}).units = ['q RAPTOR_' gdat_data.gdat_params.source{i}(1:3) ' various rho'];
gdat_data.(gdat_data.gdat_params.source{i}).dim = abcd.dim;
gdat_data.(gdat_data.gdat_params.source{i}).dimunits = abcd.dimunits;
gdat_data.(gdat_data.gdat_params.source{i}).dimunits{1} = 'various rhotor';
gdat_data.(gdat_data.gdat_params.source{i}).data_fullpath = abcd.data_fullpath;
gdat_data.(gdat_data.gdat_params.source{i}).label = ['q RAPTOR_' gdat_data.gdat_params.source{i}(1:3) ' various rho'];
end
end
end
%
% promote somethin at top level of gdat_data?
%
gdat_data.data_fullpath = {'RAP','signals_o in observer _p in predictive, one signal(q) copied at top',gdat_data.gdat_params.exp_name};
% promote predictive promoted signal if non-empty, otherwise observer one
if ~isempty(gdat_data.predictive.data)
source_promoted = 'predictive';
else
source_promoted = 'observer';
end
gdat_data.data = gdat_data.(source_promoted).data;
gdat_data.x = gdat_data.(source_promoted).x;
gdat_data.t = gdat_data.(source_promoted).t;
gdat_data.units = gdat_data.(source_promoted).units;
gdat_data.dim = gdat_data.(source_promoted).dim;
gdat_data.dimunits = gdat_data.(source_promoted).dimunits;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'rhotor', 'rhotor_edge', 'rhotor_norm', 'rhovol', 'volume_rho'}
if strcmp(upper(gdat_data.gdat_params.equil),'FPG')
gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev;
gdat_params_eff.data_request = [{'FPG'},{'Vol'},{'AUGD'}];
gdat_data = gdat_aug(gdat_data.shot,gdat_params_eff);
gdat_data.gdat_params = gdat_params_prev;
gdat_data.label = 'Vol';
gdat_data.data_fullpath = [gdat_data.data_fullpath ' yields only edge Volume from FPG/Vol'];
return
end
[gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data);
% since Lpf depends on time, need to load all first and then loop over time for easier mapping
[phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
ndimrho = size(phi_tree.data,2);
if ndimrho==NTIME_Lpf
% data seems to be transposed
ndimrho = size(phi_tree.data,1);
itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
else
itotransposeback = 0;
end
phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
ijnan=find(isnan(phi_tree.value));
phi_tree.value(ijnan)=0;
[Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
[psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
%
switch data_request_eff
case {'volume_rho', 'rhovol'}
for it=1:NTIME
Lpf1 = Lpf1_t(it);
psi_it=psi_tree.value(it,Lpf1:-1:1)';
psi_axis(it)= psi_it(1);
psi_lcfs(it)= psi_it(end);
gdat_data.x(:,it) = sqrt(abs((psi_it-psi_axis(it)) ./(psi_lcfs(it)-psi_axis(it))));
gdat_data.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)';
% gdat_data.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
gdat_data.rhovolnorm(:,it) = sqrt(abs(gdat_data.vol(:,it) ./ gdat_data.vol(end,it)));
end
gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits{1} = 'rhopolnorm';
gdat_data.dimunits{2} = 'time [s]';
gdat_data.volume_edge = gdat_data.vol(end,:);
if strcmp(data_request_eff,'volume_rho')
gdat_data.data = gdat_data.vol;
gdat_data.units = Vol.units;
else strcmp(data_request_eff,'rhovol');
gdat_data.data = gdat_data.rhovolnorm;
gdat_data.units = ' ';
end
case {'rhotor', 'rhotor_edge', 'rhotor_norm'}
b0=gdat_aug(shot,'b0');
gdat_data.b0 = interpos(b0.t,b0.data,gdat_data.t,-1);
for it=1:NTIME
Lpf1 = Lpf1_t(it);
gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)';
gdat_data.rhotornorm(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.phi(end,it)));
gdat_data.rhotor(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.b0(it) ./ pi));
psi_it=psi_tree.value(it,Lpf1:-1:1)';
psi_axis(it)= psi_it(1);
psi_lcfs(it)= psi_it(end);
gdat_data.x(:,it) = sqrt(abs((psi_it-psi_axis(it)) ./(psi_lcfs(it)-psi_axis(it))));
end
% always provide rhotor_edge so field always exists
gdat_data.rhotor_edge = gdat_data.rhotor(end,:);
if strcmp(data_request_eff,'rhotor')
gdat_data.data = gdat_data.rhotor;
gdat_data.units = 'm';
gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits{1} = 'rhopolnorm';
gdat_data.dimunits{2} = 'time [s]';
elseif strcmp(data_request_eff,'rhotor_edge')
gdat_data.data = gdat_data.rhotor(end,:);
gdat_data.units = 'm';
gdat_data.dim{1} = gdat_data.t;
gdat_data.dimunits{1} = 'time [s]';
elseif strcmp(data_request_eff,'rhotor_norm')
gdat_data.data = gdat_data.rhotornorm;
gdat_data.units = ' ';
gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits{1} = 'rhopolnorm';
gdat_data.dimunits{2} = 'time [s]';
end
end
gdat_data.data_fullpath = [DIAG '/PFL extract in .data: ' data_request_eff];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'sxr'}
% sxr from B by default or else if 'source','else' is provided
if ~isfield(gdat_data.gdat_params,'freq')|| isempty(gdat_data.gdat_params.freq)
gdat_data.gdat_params.freq = 'slow';
end
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = 'G';
elseif length(gdat_data.gdat_params.source)>1 || ...
upper(gdat_data.gdat_params.source)<'F' || upper(gdat_data.gdat_params.source)>'M'
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
dim1_len = 'from_chord_nb'; % thus up to max(chords_ok_nb)
if ~isfield(gdat_data.gdat_params,'camera') || isempty(gdat_data.gdat_params.camera)
gdat_data.gdat_params.camera = [];
elseif ischar(gdat_data.gdat_params.camera) && strcmp(gdat_data.gdat_params.camera,'central')
gdat_data.gdat_params.source = 'J';
gdat_data.gdat_params.camera = 49;
elseif isnumeric(gdat_data.gdat_params.camera)
% ok keep the array, but first dim to contain just the related chords
dim1_len='nb_of_chords';
else
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
if length(gdat_data.gdat_params.camera)==1; dim1_len='nb_of_chords'; end
gdat_data.gdat_params.source = upper(gdat_data.gdat_params.source);
%
if ~isfield(gdat_data.gdat_params,'time_interval')
gdat_data.gdat_params.time_interval = [];
end
[aa,bb]=unix(['ssh ' '-o "StrictHostKeyChecking no" sxaug20.aug.ipp.mpg.de WhichSX ' num2str(shot) ' ' ...
upper(gdat_data.gdat_params.source)]);
chords_ok_diags=regexpi(bb,'(?<chord>[A-Z]_[0-9]+)\saddress:\s[0-9]+\sDiag:\s(?<diag>[A-Z]+)','names');
if isempty(chords_ok_diags)
chords_ok_nb = [];
else
for i=1:length(chords_ok_diags)
chords_ok_nb(i) = str2num(chords_ok_diags(i).chord(3:end));
end
end
if isempty(gdat_data.gdat_params.camera);
gdat_data.gdat_params.camera = chords_ok_nb;
else
for i=1:length(gdat_data.gdat_params.camera)
ij=find(chords_ok_nb==gdat_data.gdat_params.camera(i));
if ~isempty(ij)
chords_ok_diags_tmp(i) = chords_ok_diags(ij);
else
% chord not in use
chords_ok_diags_tmp(i).chord = [];
chords_ok_diags_tmp(i).diag = [];
end
end
chords_ok_diags = chords_ok_diags_tmp;
chords_ok_nb = gdat_data.gdat_params.camera;
end
exp_name_eff = 'AUGD';
icount = 0;
nnth = 1;
if isnumeric(gdat_data.gdat_params.freq) && gdat_data.gdat_params.freq>1;
nnth = floor(gdat_data.gdat_params.freq+0.5);
gdat_data.gdat_params.freq = nnth;
end
for i=1:length(gdat_data.gdat_params.camera)
if ischar(gdat_data.gdat_params.freq) && strcmp(gdat_data.gdat_params.freq,'slow'); chords_ok_diags(i).diag = 'SSX'; end
if ~isempty(chords_ok_diags(i).diag) && ~isempty(chords_ok_diags(i).chord)
[a,e]=rdaAUG_eff(shot,chords_ok_diags(i).diag,chords_ok_diags(i).chord,exp_name_eff, ...
gdat_data.gdat_params.time_interval,gdat_data.gdat_params.extra_arg_sf2sig);
else
a.data = [];
a.t = [];
end
if ~isempty(a.data)
icount = icount + 1;
if icount == 1
% first time has data
% assume all chords have same time base even if from different shotfile
% time missing one point
if length(a.value) == length(a.t)+1
a.t=linspace(a.range(1),a.range(2),length(a.value));
a_time.index(2) = length(a.value);
end
gdat_data.t = a.t(1:nnth:end);
gdat_data.units = a.units;
if strcmp(dim1_len,'from_chord_nb')
gdat_data.data = NaN*ones(max(chords_ok_nb),length(gdat_data.t));
gdat_data.dim{1} = [1:max(chords_ok_nb)]; % simpler to have index corresponding to chord number, except for central
else
gdat_data.data = NaN*ones(length(chords_ok_nb),length(gdat_data.t));
gdat_data.dim{1} = chords_ok_nb;
end
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits = [{'chord nb'}; {'s'}];
gdat_data.data_fullpath = ['sxr from source=''' gdat_data.gdat_params.source ' uses shotfiles SX' ...
setdiff(unique(strvcat(chords_ok_diags(:).diag)),'SX')'];
gdat_data.label = ['SXR/' upper(gdat_data.gdat_params.source)];
end
try
if strcmp(dim1_len,'from_chord_nb')
gdat_data.data(chords_ok_nb(i),:) = a.data(1:nnth:end);
else
gdat_data.data(i,:) = a.data(1:nnth:end);
end
catch
if (gdat_params.nverbose>=1); disp(['problem with ' chords_ok_diags(i).diag ' ; ' chords_ok_diags(i).chord]); end
end
else
% add fields not yet defined in case all cases have empty data
end
end
gdat_data.chords = chords_ok_nb;
if length(gdat_data.dim)>=1
gdat_data.x = gdat_data.dim{1};
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'transp'}
extra_arg_sf2sig_eff_string = '';
if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
extra_arg_sf2sig_eff_string = [',' gdat_data.gdat_params.extra_arg_sf2sig];
end
% most of the times the exp for the shotfile should be provided
shotfile_exp_eff = gdat_params.exp_name;
diagname='TRA';
gdat_params.equil = 'TRE';
gdat_data.gdat_params = gdat_params;
TRANSP_signals;
for i=1:size(transp_sig,1)
if strcmp(lower(transp_sig{i,2}),'signal') || strcmp(lower(transp_sig{i,2}),'signal-group')
try
abcd = rdaAUG_eff(shot,diagname,transp_sig{i,1},shotfile_exp_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
eval(['gdat_data.' transp_sig{i,1} '= abcd;']);
% eval(['[gdat_data.' transp_sig{i,1} ',e]=rdaAUG_eff(shot,diagname,''' transp_sig{i,1} ''',shotfile_exp_eff);']);
catch
eval(['gdat_data.' transp_sig{i,1} '=[];']);
end
elseif strcmp(lower(transp_sig{i,2}),'area-base')
clear adata_area
try
eval(['[adata_area]=sf2ab(diagname,shot,transp_sig{i,1},''-exp'',shotfile_exp_eff' extra_arg_sf2sig_eff_string ');']);
adata_area.data = adata_area.value{1};
catch
adata_area.value = cell(0);
end
eval(['gdat_data.' transp_sig{i,1} '=adata_area;']);
elseif strcmp(lower(transp_sig{i,2}),'time-base')
clear adata_time
try
eval(['[adata_time]=sf2tb(diagname,shot,transp_sig{i,1},''-exp'',shotfile_exp_eff,shotfile_exp_eff' extra_arg_sf2sig_eff_string ');']);
adata_time.data = adata_time.value;
catch
adata_time.value = [];
end
eval(['gdat_data.' transp_sig{i,1} '=adata_time;']);
end
end
% copy TIME to .t
if isfield(gdat_data,'TIME') && isfield(gdat_data.TIME,'value') && ~isempty(gdat_data.TIME.value)
gdat_data.t = gdat_data.TIME.value;
gdat_data.dim{1} = gdat_data.t;
gdat_data.dimunits{1} = gdat_data.TIME.unit;
end
% add equilibrium
params_eff = gdat_params;
params_eff.data_request = 'equil';
gdat_data.equil = gdat_aug(shot,params_eff);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'zeff', 'z_eff'}
zeff_max_in_main = 10.;
zeff_min_in_main = 0.0;
sources_avail = {'ascii','idz'};
fields_to_copy = {'data','units','dim','dimunits','t','x','data_fullpath','label','help'};
for i=1:length(sources_avail)
for j=1:length(fields_to_copy)
gdat_data.(sources_avail{i}).(fields_to_copy{j}) = [];
end
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,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
%
% start with ascii then IDZ, each time copy to main, so if IDZ is present it supersedes
%
if any(strmatch('ascii',gdat_data.gdat_params.source))
% try ascii file from Rachael Mc Dermot in folder /afs/ipp-garching.mpg.de/home/r/rmcdermo/idl/Zeff/ascii
% if not present ask Rachael
dir_ascii = '/afs/ipp-garching.mpg.de/home/r/rmcdermo/idl/Zeff/ascii';
[a,b] = unix(['ls ' dir_ascii '/*' num2str(gdat_data.shot) '*.txt']);
if a==0
ij1=findstr(b,'/afs');
ij2=findstr(b,'.txt');
file_zeff = b(ij1:ij2+3);
aaa=textread(file_zeff,'%f');
if mod(numel(aaa),2)~=0
error(['problems with odd nb of points in file: ' file_zeff]);
else
nbpoints=length(aaa)/2;
end
zeff.t = aaa(1:nbpoints);
zeff.data = aaa(nbpoints+1:end);
gdat_data.ascii.data = min(max(zeff.data,zeff_min_in_main),zeff_max_in_main);
gdat_data.ascii.units = '';
gdat_data.ascii.t = zeff.t;
gdat_data.ascii.dim{1} = gdat_data.ascii.t;
gdat_data.ascii.dimunits{1} = 'time [s]';
gdat_data.ascii.data_fullpath = [' from file: ' file_zeff];
[a1,a2,a3]=fileparts(file_zeff);
gdat_data.ascii.label = [a2 a3];
for i=1:length(fields_to_copy)
if isfield(gdat_data.ascii,fields_to_copy{i})
gdat_data.(fields_to_copy{i}) = gdat_data.ascii.(fields_to_copy{i});
end
end
else
if gdat_data.gdat_params.nverbose >= 3
disp(['No ascii file for zeff in folder: ' dir_ascii '; ask Rachael if need be'])
end
end
end
%
% IDZ
%
if any(strmatch('idz',gdat_data.gdat_params.source))
gdat_data.fit.full_path = 'PROFILES from IDZ: Zeff, unc, lower, upper';
params_eff = gdat_params;
params_eff.data_request={'IDZ','Zeff', params_eff.exp_name};
try
[gdat_data_idz_zeff,params_kin,error_status]=gdat_aug(shot,params_eff);
if ~isempty(gdat_data_idz_zeff.data)
gdat_data.idz.data = gdat_data_idz_zeff.data;
gdat_data.idz.units = '';
gdat_data.idz.t = gdat_data_idz_zeff.t;
gdat_data.idz.x = gdat_data_idz_zeff.x;
gdat_data.idz.dim = gdat_data_idz_zeff.dim;
gdat_data.idz.dimunits = gdat_data_idz_zeff.dimunits;
gdat_data.idz.dimunits{1} = 'rhopol';
gdat_data.idz.data_fullpath = gdat_data_idz_zeff.data_fullpath;
gdat_data.idz.label = sprintf('%s/',gdat_data_idz_zeff.data_fullpath{:});
gdat_data.idz.label = ['mean(' gdat_data.idz.label(1:end-1) ',1)'];
fields_to_not_copy = {'data','dim','x','dimunits'}; % only cst value at this stage to main data, should heva zeff_rho to have profile
for i=1:length(fields_to_copy)
if isfield(gdat_data.idz,fields_to_copy{i}) && ~any(strmatch(fields_to_copy{i},fields_to_not_copy))
gdat_data.(fields_to_copy{i}) = gdat_data.idz.(fields_to_copy{i});
end
end
gdat_data.data = min(max(mean(gdat_data_idz_zeff.data,1)',zeff_min_in_main),zeff_max_in_main);
gdat_data.t = gdat_data_idz_zeff.t';
gdat_data.dim{1} = gdat_data.t;
gdat_data.dimunits{1} = gdat_data_idz_zeff.dimunits{2};
else
if gdat_data.gdat_params.nverbose >= 3
disp(['No data in shotfile IDZ ; ask Rainer if need be'])
end
end
catch ME
if gdat_data.gdat_params.nverbose >= 3
disp(['No shot file IDZ ; ask Rainer if need be'])
end
end
end
%
otherwise
if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_aug']); end
error_status=901;
return
end
else
if (gdat_params.nverbose>=1); warning(['AUG method=' mapping_for_AUG.method ' not known yet, contact Olivier.Sauter@epfl.ch']); end
error_status=602;
return
end
gdat_data.gdat_params.help = aug_help_parameters(fieldnames(gdat_data.gdat_params));
gdat_data.mapping_for.aug = mapping_for_aug;
gdat_params = gdat_data.gdat_params;
error_status=0;
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% functions for portions called several times
function [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data);
% get basic params to be able to read results in EQ-like shotfiles
% M_Rmesh,N_Zmesh only needed for equil when 2D quantities are required
%
exp_name_eff = gdat_data.gdat_params.exp_name;
DIAG = [];
NTIME_Lpf = [];
NTIME = [];
Lpf1_t = [];
Lpf_SOL = [];
M_Rmesh = [];
N_Zmesh = [];
extra_arg_sf2sig_eff_string = '';
if ~strcmp(gdat_data.gdat_params.extra_arg_sf2sig,'[]')
extra_arg_sf2sig_eff_string = [',' gdat_data.gdat_params.extra_arg_sf2sig];
end
shot=gdat_data.shot;
exp_name_eff = gdat_data.gdat_params.exp_name;
gdat_data.gdat_params.equil = upper(gdat_data.gdat_params.equil);
DIAG = gdat_data.gdat_params.equil; % 'EQI' by default at this stage, should become EQH?
% eval(['M_Rmesh_par = sf2par(DIAG,shot,''M'',''PARMV''' extra_arg_sf2sig_eff_string ');']);
[M_Rmesh_par,e]=rdaAUG_eff(shot,DIAG,'PARMV',gdat_data.gdat_params.exp_name,[], ...
extra_arg_sf2sig_eff_string,['param:' 'M']);
M_Rmesh = M_Rmesh_par.value + 1; % nb of points
% eval(['N_Zmesh_par = sf2par(DIAG,shot,''N'',''PARMV''' extra_arg_sf2sig_eff_string ');']);
[N_Zmesh_par,e]=rdaAUG_eff(shot,DIAG,'PARMV',gdat_data.gdat_params.exp_name,[], ...
extra_arg_sf2sig_eff_string,['param:' 'N']);
N_Zmesh = N_Zmesh_par.value + 1; % nb of points
% eval(['NTIME_par = sf2par(DIAG,shot,''NTIME'',''PARMV''' extra_arg_sf2sig_eff_string ');']);
[NTIME_par,e]=rdaAUG_eff(shot,DIAG,'PARMV',gdat_data.gdat_params.exp_name,[], ...
extra_arg_sf2sig_eff_string,['param:' 'NTIME']);
NTIME = NTIME_par.value; % nb of points
Lpf_par = rdaAUG_eff(shot,DIAG,'Lpf',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
% since June, nb of time points in EQ results is not consistent with NTIME and time
% It seems the first NTIME points are correct, so use this explicitely
NTIME_Lpf = length(Lpf_par.value);
if ~isempty(NTIME) && ~isempty(NTIME_Lpf) && isnumeric(NTIME) && isnumeric(NTIME_Lpf) && (NTIME < NTIME_Lpf)
if gdat_data.gdat_params.nverbose>=3; disp('WARNING: nb of times points smaller then equil results, use first NTIME points'); end
elseif ischar(NTIME) || ischar(NTIME_Lpf) || (NTIME > NTIME_Lpf)
if gdat_data.gdat_params.nverbose >= 1
if ischar(NTIME) || ischar(NTIME_Lpf)
disp(['probably no data, NTIME_Lpf = ' NTIME_Lpf])
else
disp('ERROR: nb of times points LARGER then equil results')
disp('this is unexpected, so stop there and ask Olivier.Sauter@epfl.ch')
end
end
return
end
Lpf_tot = Lpf_par.value(1:NTIME); % nb of points: 100000*nb_SOL points + nb_core
Lpf_SOL = fix(Lpf_tot/100000);
Lpf1_t = mod(Lpf_tot,100000)+1; % nb of Lpf points
[equil_time,e]=rdaAUG_eff(shot,DIAG,'time',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
gdat_data.t = equil_time.value(1:NTIME);