From a5ea0a157c6c3b1eb5264323ce4a6d601288ba53 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Wed, 15 Jan 2020 16:49:56 +0100
Subject: [PATCH] add machine='imas' and getting 'ids','source','ids_name' from
IMAS database
---
matlab/IMAS/gdat_imas.m | 716 ++++++++++++++++++++++++++++
matlab/IMAS/imas_help_parameters.m | 60 +++
matlab/IMAS/imas_requests_mapping.m | 43 ++
3 files changed, 819 insertions(+)
create mode 100644 matlab/IMAS/gdat_imas.m
create mode 100644 matlab/IMAS/imas_help_parameters.m
create mode 100644 matlab/IMAS/imas_requests_mapping.m
diff --git a/matlab/IMAS/gdat_imas.m b/matlab/IMAS/gdat_imas.m
new file mode 100644
index 00000000..1547090f
--- /dev/null
+++ b/matlab/IMAS/gdat_imas.m
@@ -0,0 +1,716 @@
+function [gdat_data,gdat_params,error_status,varargout] = gdat_imas(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
+%
+% IMAS means assume access data from IDS database ala IMAS with tokamakname and dataversion and run_number using imas_open_env and ids_get calls
+% Should be the same on ITER hpc and gateway, to be checked
+%
+% At first only "ids" request implemented, but keywords like ip etc should relate to specific IDS as for any machine
+%
+% Inputs:
+%
+% no inputs: return default parameters in a structure form in gdat_params
+% shot: shot number (if empty uses 'idx_imas_open_env' in params from previous call or input param)
+% 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','IMAS' (the default machine is the local machine)
+%
+%
+% Outputs:
+%
+% gdat_data: structure containing the data, the time trace if known and other useful information
+% gdat_data.t : time trace
+% gdat_data.data: requested data values
+% gdat_data.dim : values of the various coordinates related to the dimensions of .data(:,:,...)
+% note that one of the dim is the time, replicated in .t for clarity
+% gdat_data.dimunits : units of the various dimensions, 'dimensionless' if dimensionless
+% gdat_data.error_bar : if provided
+% gdat_data.gdat_call : list of parameters provided in the gdat call (so can be reproduced)
+% gdat_data.shot: shot number
+% gdat_data.machine: machine providing the data
+% gdat_data.gdat_request: keyword for gdat if relevant
+% gdat_data.data_fullpath: full path to the data node if known and relevant, or relevant expression called if relevant
+% gdat_data.gdat_params: copy gdat_params for completeness (gdat_params contains a .help structure detailing each parameter)
+% gdat_data.xxx: any other relevant information
+%
+% gdat_params contains the options relevant for the called data_request. It also contains a help structure for each option
+% eg.: param1 in gdat_params.param1 and help in gdat_params.help.param1
+%
+% Note for liuqe parameter: At this time we have liuqe_fortran (which used to be the default) and liuqe_matlab (which is new and now becomes/is the default)
+% We still have liuqe1, 2 and 3 nodes for either of them, so you can choose:
+% 'liuqe',1 (2 or 3): to get the default liuqe 1, 2 or 3 (which is now the matlab version)
+% 'liuqe',11 (12 or 13): to get liuqe_fortran nodes 1, 2 or 3
+% 'liuqe',21 (22 or 23): to get liuqe_matlab nodes 1, 2 or 3 (may be needed if we set the default to liuqe_fortran for old shots)
+% 'liuqe',-1 : to get FBTE result
+%
+% Examples:
+% (should add working examples for various machines (provides working shot numbers for each machine...))
+%
+% to get the list of ids: aa=gdat([],'ids','machine','imas'); aa.gdat_params.sources_available
+% to get an empty ids: aa=gdat([],'ids','machine','imas','source','core_profiles'); aa.core_profiles
+%
+% [a1,a2]=gdat;
+% a2.data_request = 'ids'; a2.source='equilibrium'; a2.machine = 'imas';
+% a3=gdat(48836,a2); % gives input parameters as a structure, allows to call the same for many shots
+% a4=gdat('opt1',123,'opt2',[1 2 3],'shot',48832,'data_request','Ip','opt3','aAdB'); % all in pairs
+% a5=gdat(48836,'ids','source','core_profiles'); % standard call for a given IDS core_profiles, use {'',''} for a list of IDSs
+% a6=gdat(48836,'ip'); % standard call for a specific keyword
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%
+% Comments for local developer:
+% gdat is just a "header routine" calling the gdat for the specific machine gdat_`machine`.m which can be called
+% directly, thus which should be able to treat the same type of input arguments
+%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Prepare some variables, etc
+
+varargout{1}=cell(1,1);
+error_status=1;
+
+% construct main default parameters structure
+gdat_params.data_request = '';
+default_machine = 'imas';
+
+gdat_params.machine=default_machine;
+gdat_params.doplot = 0;
+gdat_params.nverbose = 1;
+gdat_params.run = []; % should be defined by user to make sure
+gdat_params.occurence = 0; % default/usual value
+
+% special case need to treat here if params given as structure in data_request:
+if isstruct(data_request)
+ gdat_params = data_request;
+end
+
+% construct list of keywords from global set of keywords and specific IMAS 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 IMAS specific to all:
+if ~isempty(data_request_names.imas)
+ imas_names = fieldnames(data_request_names.imas);
+ for i=1:length(imas_names)
+ data_request_names.all.(imas_names{i}) = data_request_names.imas.(imas_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
+ if ischar(data_request);
+ data_request_eff = data_request;
+ data_request = data_request; % should not lower until see if keyword
+ elseif isstruct(data_request)
+ data_request_eff = data_request.data_request;
+ end
+end
+gdat_data.gdat_request = data_request_names_all; % so if return early gives list of possible request names
+gdat_data.gdat_params.help = imas_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)
+ % if shot already open, will use idx from gdat_params, otherwise returns
+ do_mdsopen_mdsclose = 0;
+ if ~isfield(gdat_data.gdat_params,'idx_imas_open_env')
+ if ~strcmp(data_request_eff,'ids'); return; end % empty shot return empty ids so valid command
+ end
+ elseif isnumeric(shot)
+ gdat_data.shot = shot;
+ elseif ischar(shot)
+ ivarargin_first_char = 1;
+ else
+ if gdat_params.nverbose>=1; warning('type of 1st argument unexpected, should be numeric or char'); end
+ error_status=2;
+ return
+ end
+ if nargin==1
+ % Only shot number given. If there is a default data_request set it and continue, otherwise return
+ return
+ end
+end
+% 2nd input argument if not part of pairs
+if nargin>=2 && ivarargin_first_char~=1
+ if isempty(data_request)
+ return
+ end
+ % 2nd arg can be a structure with all options except shot_number, or a string for the pathname or keyword, or the start of pairs string/value for the parameters
+ if isstruct(data_request)
+ if ~isfield(data_request,'data_request')
+ if gdat_params.nverbose>=1; warning('expects field data_request in input parameters structure'); end
+ error_status=3;
+ return
+ end
+ data_request_eff = data_request.data_request;
+ % some request are case sensitive because mds names can be like \magnetics::ipol[*,"OH_001"]
+ data_request.data_request = data_request.data_request;
+ gdat_params = data_request;
+ else
+ % since data_request is char check from nb of args if it is data_request or a start of pairs
+ if mod(nargin-1,2)==0
+ ivarargin_first_char = 2;
+ else
+ ivarargin_first_char = 3;
+ data_request_eff = data_request;
+ end
+ end
+end
+
+if ~isstruct(data_request)
+ gdat_params.data_request = data_request_eff;
+end
+
+% if start pairs from shot or data_request, shift varargin
+if ivarargin_first_char==1
+ varargin_eff{1} = shot;
+ varargin_eff{2} = data_request;
+ varargin_eff(3:nargin) = varargin(:);
+elseif ivarargin_first_char==2
+ varargin_eff{1} = data_request;
+ varargin_eff(2:nargin-1) = varargin(:);
+else
+ varargin_eff(1:nargin-2) = varargin(:);
+end
+
+% extract parameters from pairs of varargin:
+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})
+ gdat_params.(lower(varargin_eff{i})) = varargin_eff{i+1}; % cannot enforce lower(params value) if filename or else are case sensitive
+ else
+ if gdat_params.nverbose>=1; warning(['input argument nb: ' num2str(i) ' is incorrect, expects a character string']); end
+ error_status=401;
+ return
+ end
+ end
+ else
+ if gdat_params.nverbose>=1; warning('number of input arguments incorrect, cannot make pairs of parameters'); end
+ error_status=402;
+ return
+ end
+end
+
+% if it is a request_keyword copy it:
+if ischar(data_request_eff) || length(data_request_eff)==1
+ ij=strmatch(lower(data_request_eff),data_request_names_all,'exact');
+else
+ ij=[];
+end
+if ~isempty(ij);
+ gdat_data.gdat_request = data_request_names_all{ij};
+ gdat_params.data_request = gdat_data.gdat_request;
+ if isfield(data_request_names.all.(data_request_names_all{ij}),'description') && ~isempty(data_request_names.all.(data_request_names_all{ij}).description)
+ % copy description of keyword
+ gdat_data.request_description = data_request_names.all.(data_request_names_all{ij}).description;
+ end
+end
+
+% special treatment if shot and data_request given within pairs
+if isfield(gdat_params,'shot')
+ shot = gdat_params.shot; % should use only gdat_params.shot but change shot to make sure
+ gdat_data.shot = gdat_params.shot;
+ gdat_params=rmfield(gdat_params,'shot');
+end
+
+if ~isfield(gdat_params,'data_request') || isempty(gdat_params.data_request)
+ % warning('input for ''data_request'' missing from input arguments') % might be correct, asking for list of requests
+ error_status=5;
+ return
+end
+gdat_data.gdat_params = gdat_params; % This means that before here it is gdat_params which should be updated
+clear gdat_params; % to make sure one uses gdat_data.gdat_params now
+
+% 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;
+if ~isfield(gdat_data.gdat_params,'database_user')
+ gdat_data.gdat_params.database_user = getenv('USER');
+ if gdat_data.gdat_params.nverbose >= 3; disp('set database_user to username as default since not provided'); end
+end
+if ~isfield(gdat_data.gdat_params,'tokamak_name')
+ gdat_data.gdat_params.tokamak_name = gdat_data.gdat_params.machine;
+ if gdat_data.gdat_params.nverbose >= 3; disp('set tokamak_name to machine as default since not provided'); end
+end
+if ~isfield(gdat_data.gdat_params,'data_major_version')
+ gdat_data.gdat_params.data_major_version = '3';
+ if gdat_data.gdat_params.nverbose >= 3; disp('set data_major_version to ''3'' as default since not provided'); end
+end
+error_status = 6; % at least reached this level
+gdat_params = gdat_data.gdat_params;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Specifications on how to get the data provided in imas_requests_mapping
+mapping_for_imas = imas_requests_mapping(data_request_eff,shot);
+gdat_data.label = mapping_for_imas.label;
+ishot=NaN;
+if do_mdsopen_mdsclose
+ if ~isfield(gdat_data.gdat_params,'run') || isempty(gdat_data.gdat_params.run)
+ warning('run number not provided, cannot continue')
+ return
+ end
+ try
+ gdat_data.gdat_params.idx_imas_open_env = imas_open_env('ids', shot, gdat_data.gdat_params.run,gdat_data.gdat_params.database_user, ...
+ gdat_data.gdat_params.tokamak_name,gdat_data.gdat_params.data_major_version);
+ disp('')
+ catch ME
+ warning('could not imas_open_env with following gdat_data.gdat_params:');
+ gdat_data.gdat_params
+ rethrow(ME)
+ end
+ if gdat_data.gdat_params.idx_imas_open_env < 0
+ if gdat_data.gdat_params.nverbose>=1
+ warning(['idx < 0, shot/run not opened properly with gdat_data.gdat_params:']);
+ gdat_data.gdat_params
+ end
+ return
+ end
+else
+ if isfield(gdat_data.gdat_params,'idx_imas_open_env')
+ if gdat_data.gdat_params.idx_imas_open_env < 0
+ if gdat_data.gdat_params.nverbose>=1
+ warning(['idx < 0, shot/run not opened properly with gdat_data.gdat_params:']);
+ gdat_data.gdat_params
+ end
+ return
+ end
+ else
+ if ~isempty(shot)
+ if gdat_data.gdat_params.nverbose>=1
+ warning(['no idx in gdat_params, cannot open shot with gdat_data.gdat_params:']);
+ gdat_data.gdat_params
+ end
+ return
+ end
+ end
+end
+gdat_params = gdat_data.gdat_params;
+
+% fill again at end to have full case, but here to have present status in case of early return
+gdat_data.gdat_params.help = imas_help_parameters(fieldnames(gdat_data.gdat_params));
+gdat_data.mapping_for.imas = mapping_for_imas;
+gdat_params = gdat_data.gdat_params;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if strcmp(mapping_for_imas.method,'expression')
+ % 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_imas.expression ';'];
+ eval([mapping_for_imas.expression ';']);
+ if isempty(gdat_tmp) || (~isstruct(gdat_tmp) & ~isobject(gdat_tmp))
+ if (gdat_params.nverbose>=1); warning(['expression does not create a gdat_tmp structure: ' mapping_for_imas.expression]); end
+ error_status=801;
+ return
+ end
+ tmp_fieldnames = setdiff(fieldnames(gdat_tmp),{'gdat_request','label'}); % could/should also remove label in any case
+ if sum(strcmp(tmp_fieldnames,'data'))==0 % note: cannot do isfield since gdat_tmp might be an object
+ if (gdat_params.nverbose>=1); warning(['expression does not return a child name ''data'' for ' data_request_eff]); end
+ end
+ for i=1:length(tmp_fieldnames)
+ gdat_data.(tmp_fieldnames{i}) = gdat_tmp.(tmp_fieldnames{i});
+ end
+ % add .t and .x in case only dim is provided
+ % do not allow shifting of timedim since should be treated in the relevant expression
+ ijdim=find(strcmp(tmp_fieldnames,'dim')==1);
+ if ~isempty(ijdim)
+ nbdims = length(gdat_data.dim);
+ if mapping_for_imas.timedim==-1;
+ mapping_for_imas.timedim = nbdims;
+ if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_imas.timedim = nbdims-1; end
+ end
+ dim_nontim = setdiff([1:nbdims],mapping_for_imas.timedim);
+ ijt=find(strcmp(tmp_fieldnames,'t')==1);
+ if isempty(ijt)
+ gdat_data.t = gdat_data.dim{mapping_for_imas.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_imas.expression;
+ if isfield(gdat_tmp,'help')
+ gdat_data.help = gdat_tmp.help;
+ else
+ gdat_data.help = [];
+ end
+ if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
+ [gdat_data] = gdat2time_out(gdat_data,1);
+ end
+ end
+ % end of method "expression"
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+elseif strcmp(mapping_for_imas.method,'switchcase')
+ switch data_request_eff % not lower(...) since data_request_eff should be lower case already at this stage
+ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % First the request names valid for "all" machines:
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'ids'}
+ ids_empty_path = fullfile(fileparts(mfilename('fullpath')),'..','IMAS_IMAS','ids_empty');
+
+ params_eff = gdat_data.gdat_params;
+ if isfield(params_eff,'source') && ~isempty(params_eff.source)
+ ids_top_name = params_eff.source;
+ else
+ ids_top_name = [];
+ warning('gdat:EmptyIDSName','Need an ids name in ''source'' parameter\n check substructure gdat_params.sources_available for an ids list');
+ addpath(ids_empty_path);
+ assert(~~exist('ids_list','file'),'could not find ids_list.m in %s',ids_empty_path);
+ gdat_data.gdat_params.sources_available = ids_list;
+ rmpath(ids_empty_path);
+ return
+ end
+ ids_gen_ok = ~~exist('ids_gen','file');
+
+ if ids_gen_ok
+ ids_empty = ids_gen(ids_top_name); % generate ids from gateway function ids_gen
+ else
+ % load empty ids structure from template file
+ fname = sprintf('ids_empty_%s',ids_top_name);
+
+ try
+ assert(~~exist(ids_empty_path,'dir'),'folder %s not found',ids_empty_path);
+ addpath(ids_empty_path);
+ assert(~~exist(fname,'file'),'file %s does not exist in %s',fname,ids_empty_path);
+
+ ids_empty = eval(fname);
+ disp(['use structure in .mat file: ' ids_empty_path '/' fname]);
+ rmpath(ids_empty_path);
+ catch ME
+ fprintf('Could not load empty template for %s\n',ids_top_name);
+ fprintf('Available templates:\n');
+ disp(dir(ids_empty_path));
+ rmpath(ids_empty_path);
+ rethrow(ME);
+ end
+ end
+ if ~isfield(gdat_data.gdat_params,'error_bar') || isempty(gdat_data.gdat_params.error_bar)
+ gdat_data.gdat_params.error_bar = 'delta';
+ end
+ if ~isfield(gdat_data.gdat_params,'cocos_in') || isempty(gdat_data.gdat_params.cocos_in)
+ gdat_data.gdat_params.cocos_in = 11;
+ end
+ if ~isfield(gdat_data.gdat_params,'cocos_out') || isempty(gdat_data.gdat_params.cocos_out)
+ gdat_data.gdat_params.cocos_out = 11;
+ end
+ if ~isfield(gdat_data.gdat_params,'ipsign_out') || isempty(gdat_data.gdat_params.ipsign_out)
+ gdat_data.gdat_params.ipsign_out = 0;
+ end
+ if ~isfield(gdat_data.gdat_params,'b0sign_out') || isempty(gdat_data.gdat_params.b0sign_out)
+ gdat_data.gdat_params.b0sign_out = 0;
+ end
+ if ~isfield(gdat_data.gdat_params,'ipsign_in') || isempty(gdat_data.gdat_params.ipsign_in)
+ if gdat_data.gdat_params.ipsign_out~=0
+ gdat_data.gdat_params.ipsign_in = mdsvalue('\pcs::data:iohfb');
+ else
+ gdat_data.gdat_params.ipsign_in = 0;
+ end
+ end
+ if ~isfield(gdat_data.gdat_params,'b0sign_in') || isempty(gdat_data.gdat_params.b0sign_in)
+ if gdat_data.gdat_params.b0sign_out~=0
+ gdat_data.gdat_params.b0sign_in = mdsvalue('\pcs::data:if36fb');
+ else
+ gdat_data.gdat_params.b0sign_in = 0;
+ end
+ end
+ try
+ if isfield(gdat_data.gdat_params,'idx_imas_open_env') && gdat_data.gdat_params.idx_imas_open_env >= 0
+ if gdat_data.gdat_params.occurence == 0
+ ids_top = ids_get(gdat_data.gdat_params.idx_imas_open_env,ids_top_name);
+ else
+ ids_top = ids_get(gdat_data.gdat_params.idx_imas_open_env,[ids_top_name '/' num2str(gdat_data.gdat_params.occurence)]);
+ end
+ gdat_data.(ids_top_name) = ids_top;
+ else
+ gdat_data.(ids_top_name) = ids_empty;
+ end
+ catch ME_imas_ids_get
+ disp(['there is a problem with: imas_get_ids_' ids_top_name ...
+ ' , may be check if it exists in your path or test it by itself'])
+ gdat_data.(ids_top_name) = ids_empty;
+ gdat_data.([ids_top_name '_help']) = getReport(ME_imas_ids_get);
+ rethrow(ME_imas_ids_get)
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ otherwise
+ if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_imas']); end
+ error_status=901;
+ return
+ end
+
+else
+ if (gdat_params.nverbose>=1); warning(['IMAS method=' mapping_for_imas.method ' not known yet, contact Olivier.Sauter@epfl.ch']); end
+ error_status=602;
+ return
+end
+
+gdat_data.gdat_params.help = imas_help_parameters(fieldnames(gdat_data.gdat_params));
+gdat_data.mapping_for.imas = mapping_for_imas;
+gdat_params = gdat_data.gdat_params;
+error_status=0;
+
+return
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [firthomratio] = get_fir_thom_rat_data(shot,maintracename,timebase,nverbose);
+%
+% since depends on shot number for using auto fit and thomson or thomson edge, use tracename and function here
+%
+% maintracename = 'thomson' or 'thomson_edge
+%
+% return fir_to_thomson ratio on time=timebase
+%
+% normally should use "main" thomson one built from auto fit if possible
+% take edge one if need be
+%
+% default: maintracename = "thomson"
+%
+maintracename_eff = 'thomson';
+if exist('maintracename') && ~isempty(maintracename)
+ maintracename_eff = maintracename;
+end
+
+firthomratio = NaN;
+if ~exist('timebase') || isempty(timebase)
+ if (nverbose>=1)
+ disp('need a timebase in get_fir_thom_rat_data')
+ end
+ return
+end
+firthomratio = NaN*ones(size(timebase));
+
+if strcmp(maintracename_eff,'thomson')
+ if shot>=23801
+ tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
+ if isempty(tracefirrat.data) || ischar(tracefirrat.data)
+ if (nverbose>=1); disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat'); end
+ tracefirrat=tdi('\results::thomson:fir_thom_rat');
+ end
+ else
+ tracefirrat=tdi('\results::thomson:fir_thom_rat');
+ if isempty(tracefirrat.data) || ischar(tracefirrat.data)
+ if (nverbose>=1); disp('problem with \results::thomson:fir_thom_rat: empty'); end
+ end
+ end
+elseif strcmp(maintracename_eff,'thomson_edge')
+ if shot>=23801
+ tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
+ if isempty(tracefirrat.data) || ischar(tracefirrat.data)
+ if (nverbose>=1); disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty, use thomson:fir_thom_rat'); end
+ tracefirrat=tdi('\results::thomson:fir_thom_rat');
+ end
+ else
+ tracefirrat=tdi('\results::thomson_edge:fir_thom_rat');
+ if isempty(tracefirrat.data) || ischar(tracefirrat.data)
+ if (nverbose>=1); disp('problem with \results::thomson_edge:fir_thom_rat: empty'); end
+ end
+ end
+else
+ if (nverbose>=1); disp('bad input in get_fir_thom_rat_data'); end
+ return
+end
+
+if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) && any(isfinite(tracefirrat.data)) ...
+ && ~isempty(tracefirrat.dim) && ~isempty(tracefirrat.dim{1})
+ firthomratio = interp1(tracefirrat.dim{1},tracefirrat.data,timebase);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,doedge,nverbose);
+%
+try
+ time=mdsdata('\results::thomson:times');
+catch
+ gdat_data.error_bar = [];
+ if strcmp(data_request_eff(1:2),'ne')
+ tracefirrat_data = [];
+ gdat_data.firrat=tracefirrat_data;
+ gdat_data.data_raw = gdat_data.data;
+ gdat_data.error_bar_raw = gdat_data.error_bar;
+ end
+
+ if (nverbose>=1) && shot<100000
+ warning('Problems with \results::thomson:times')
+ warning(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+ end
+ return
+end
+if isempty(time) || ischar(time)
+ thomsontimes=time;
+ if (nverbose>=1) && shot<100000
+ warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check')
+ disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+ end
+ return
+end
+edge_str_ = '';
+edge_str_dot = '';
+if doedge
+ edge_str_ = '_edge';
+ edge_str_dot = '.edge';
+end
+
+nodenameeff = ['\results::thomson' edge_str_dot ':' data_request_eff(1:2)];
+tracetdi=tdi(nodenameeff);
+if isempty(tracetdi.data) || ischar(tracetdi.data) || isempty(tracetdi.dim)
+ gdat_data.error_bar = [];
+ gdat_data.firrat = [];
+ gdat_data.data_raw = [];
+ gdat_data.error_bar_raw = [];
+ return
+end
+
+gdat_data.data=tracetdi.data'; % Thomson data as (t,z)
+tracestd=tdi(['\results::thomson' edge_str_dot ':' data_request_eff(1:2) ':error_bar']);
+gdat_data.error_bar=tracestd.data';
+gdat_data.data_fullpath=[nodenameeff '; error_bar'];
+% add fir if ne requested
+if strcmp(data_request_eff(1:2),'ne')
+ tracefirrat_data = get_fir_thom_rat_data(shot,['thomson' edge_str_],time,nverbose);
+ gdat_data.firrat=tracefirrat_data;
+ gdat_data.data_raw = gdat_data.data;
+ gdat_data.data = gdat_data.data_raw * diag(tracefirrat_data);
+ gdat_data.error_bar_raw = gdat_data.error_bar;
+ gdat_data.error_bar = gdat_data.error_bar_raw * diag(tracefirrat_data);
+ gdat_data.data_fullpath=[gdat_data.data_fullpath '; fir_thom_rat ; _raw without firrat'];
+ ij=find(isfinite(tracefirrat_data));
+ if isempty(ij)
+ gdat_data.data_fullpath=[gdat_data.data_fullpath ' WARNING use ne_raw since firrat has NaNs thus ne=NaNs; fir_thom_rat ; _raw without firrat'];
+ disp('***********************************************************************')
+ disp('WARNING: ne from Thomson has fir_thom_rat with Nans so only ne_raw can be used');
+ disp('***********************************************************************')
+ end
+end
+z=mdsdata(['\diagz::thomson_set_up' edge_str_dot ':vertical_pos']);
+gdat_data.dim=[{z};{time}];
+gdat_data.dimunits=[{'Z [m]'} ; {'time [s]'}];
+gdat_data.x=z;
+gdat_data.t=time;
+% isfield does not work since tracetdi is not a 'struct' but a tdi object, thus use fieldnames
+if any(strcmp(fieldnames(tracetdi),'units'))
+ gdat_data.units=tracetdi.units;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, system_str_dot, psitbx_str, gdat_params )
+% Input arguments:
+% gdat_data: current structure containing TS data
+% zshift: vertical shift of the equilibrium, can be a vector for time-dependent data
+% system_str_dot: '' or '.edge' indicating thomson syste,
+% psitbx_str: source for the equilibrium map, psitbx-style
+% gdat_params: parameters for the gdat call % TODO: Why is it different from gdat_data.gdat_params?
+
+% Use psitbx
+t_th = gdat_data.t;
+th_z = mdsdata(['dim_of(\thomson' system_str_dot ':te,1)']); % TODO: Isn't this gdat_data.dim{1}?
+th_r = 0.9*ones(size(th_z));
+ntt = numel(t_th);
+nch = numel(th_z); % number of chords
+
+[t_psi,status] = mdsvalue('dim_of(imas_eq("i_p",$))',psitbx_str); % TODO: replace with time_psi when imas_eq supports it
+if ~rem(status,2)
+ warning('problem with getting equilibrium time base in gdat_imas')
+ display(t_psi);
+ return
+end
+
+zshifteff=zshift;
+% set zshifteff time array same as t_psi
+switch numel(zshifteff)
+ case 1
+ zshifteff=zshifteff * ones(size(t_th));
+ case numel(t_psi)
+ zshifteff=interp1(t_psi,zshifteff,t_th);
+ case numel(t_th)
+ % ok
+ otherwise
+ if (gdat_params.nverbose>=1)
+ disp(' bad time dimension for zshift')
+ disp(['it should be 1 or ' num2str(numel(t_th)) ' or ' num2str(numel(t_psi))])
+ return
+ end
+end
+
+% Get LIUQE times corresponding to TS times
+t_eq = t_psi(iround(t_psi,t_th));
+
+% Get PSI map
+psi = psitbximas(gdat_data.shot,t_eq,'*0',psitbx_str);
+% PSITBXIMAS will have removed duplicate times
+i_psi = iround(psi.t,t_eq); % Mapping from psi times to TS times
+
+% Define grid points where to evaluate PSI
+if ~any(diff(zshifteff))
+ % Constant vector
+ pgrid = psitbxgrid('cylindrical','points',{th_r,th_z - zshifteff(1),0});
+else
+ % We need as many time points in psi as TS times
+ psi = subs_time(psi,i_psi); % Now has ntt times
+ i_psi = 1:ntt; % Trivial mapping from psi times to TS times
+ % Time-varying vector
+ th_z_eff = repmat(th_z(:),1,ntt) - repmat(reshape(zshifteff,1,ntt),nch,1);
+ th_r_eff = repmat(th_r(:),1,ntt);
+ pgrid = psitbxgrid('cylindrical','time-points',{th_r_eff,th_z_eff,0});
+end
+
+% Compute psi at TS positions
+psi_ts = psitbxf2f(psi,pgrid);
+psiscatvol.data = squeeze(psi_ts.x(:,i_psi));
+psiscatvol.dim{1} = gdat_data.x;
+psiscatvol.dim{2} = t_th;
+% NOTE: we should probably not include time points where equilibrium time is far from TS time.
+
+% Compute psi_axis at TS times
+% Do not use psi_axis node because:
+% - For legacy LIUQE, psi_axis is not computed by psitbx, so we could get a
+% different answer and cause complex numbers computing rho_psi
+% - For MATLAB LIUQE, psi_axis is only the result of asxymex whereas the
+% psitbx adds some Newton iterations so again complex numbers are possible
+psi_norm = psitbxp2p(psi,'01');
+psi_max.data = psi_norm.psimag(i_psi);
+psi_max.dim = {t_th};
diff --git a/matlab/IMAS/imas_help_parameters.m b/matlab/IMAS/imas_help_parameters.m
new file mode 100644
index 00000000..bbf4bbd7
--- /dev/null
+++ b/matlab/IMAS/imas_help_parameters.m
@@ -0,0 +1,60 @@
+function help_struct = imas_help_parameters(parameter_list)
+%
+% retrieve from present table the relevant help lines for the parameter_list{:}
+% should come from sqlite database at some point...
+%
+% return the whole help structure if parameter_list empty or not provided
+%
+% do:
+% help_struct = imas_help_parameters(fieldnames(gdat_data.gdat_params));
+%
+% to get relevant help description
+%
+
+% Defaults
+help_struct_all = struct(...
+ 'data_request', ['automatically filled in by gdat, name of request used in gdat call.' char(10) ...
+ 'contains current list of keywords if gdat called with no arguments: aa=gdat;' char(10) ...
+ 'Note shot value should not be in params so params can be used to load same data from another shot'] ...
+ ,'machine', 'machine name like ''TCV'', ''AUG'', ''IMAS'' case insensitive' ...
+ ,'doplot', '0 (default), if 1 calls gdat_plot for a new figure, -1 plot over current figure with hold all, see gdat_plot for details' ...
+ ,'nverbose','1 (default) displays warnings, 0: only errors, >=3: displays all extra information' ...
+ );
+
+% IMAS related
+help_struct_all.time = 'eqdsk: time(s) value(s) requested, by default time=1.0s (see time_out for other requests)';
+help_struct_all.time_out = ['requested time for output: data points within interval if time_out=[t1 t2], otherwise assumes series of points, uses linear interpolation in that case (default [-Inf Inf])'...
+ char(10) 'for sxr, mpx: only time interval provided in time_out is relevant'];
+help_struct_all.cocos_in_out = ['cocos_in or cocos_out value input or desired in output, uses eqdsk_cocos_transform. Note should use latter if a specific Ip and/or B0 sign' ...
+ 'is wanted. See O. Sauter et al Comput. Phys. Commun. 184 (2013) 293'];
+help_struct_all.source = sprintf('%s\n',...
+ 'ids_names for request ''ids'' like magnetics, equilibrium, etc');
+help_struct_all.error_bar = sprintf('%s\n','for ids: choice of nodes fill in and how:', ...
+ '''delta'' (default): only upper fill in with the abs(value) to add or subtract to data to get upper and lower values (symmetric)', ...
+ '''delta_with_lower'': same as delta but fill in lower node as well (with delta as well, same as upper)', ...
+ '''added'': add the delta values (old cpo style), so upper=data+error_bar and lower=data+error_bar');
+help_struct_all.write = 'eqdsk: write eqdsk while loading data (1, default) or not (0)';
+%help_struct_all. = '';
+
+if ~exist('parameter_list') || isempty(parameter_list)
+ help_struct = help_struct_all;
+ return
+end
+
+if iscell(parameter_list)
+ parameter_list_eff = parameter_list;
+elseif ischar(parameter_list)
+ % assume only one parameter provided as char string
+ parameter_list_eff{1} = parameter_list;
+else
+ disp(['unknown type of parameter_list in imas_help_parameters.m'])
+ parameter_list
+ return
+end
+
+fieldnames_all = fieldnames(help_struct_all);
+for i=1:length(parameter_list_eff)
+ if any(strcmp(fieldnames_all,parameter_list_eff{i}))
+ help_struct.(parameter_list_eff{i}) = help_struct_all.(parameter_list_eff{i});
+ end
+end
diff --git a/matlab/IMAS/imas_requests_mapping.m b/matlab/IMAS/imas_requests_mapping.m
new file mode 100644
index 00000000..575e4649
--- /dev/null
+++ b/matlab/IMAS/imas_requests_mapping.m
@@ -0,0 +1,43 @@
+function mapping = imas_requests_mapping(data_request,shot)
+%
+% Information pre-defined for gdat_imas, you can add cases here to match official cases in gdat_imas, allowing backward compatibility
+%
+% give the shot number in case data origin depends on the shot number, allows to adapt easily
+%
+
+% Defaults
+mapping = struct(...
+ 'label', '', ...
+ 'method', '', ...
+ 'expression','', ...
+ 'timedim', -1, ... % dim which is the time is the database, to copy in .t, the other dims are in .x (-1 means last dimension)
+ 'gdat_timedim',[], ... % if need to reshape data and dim orders to have timedim as gdat_timedim (shifting time to gdat_timedim)
+ 'min', -inf, ...
+ 'max', inf);
+% Note that gdat_timedim is set to timedim at end of this function if empty
+% gdat_timedim should have effective index of the time dimension in gdat
+
+if ~exist('data_request') || isempty(data_request)
+ return
+end
+
+% default label: data_request keyword itself
+mapping.label = data_request;
+
+%
+% label is used for plotting
+switch lower(data_request)
+ case 'ids'
+ mapping.timedim = 1;
+ mapping.label = 'ids ala imas';
+ mapping.method = 'switchcase';
+ otherwise
+ mapping.label = data_request;
+ mapping.method = 'switchcase';
+ mapping.expression = data_request;
+
+end
+
+if isempty(mapping.gdat_timedim)
+ mapping.gdat_timedim = mapping.timedim;
+end
--
GitLab