From 3b61bb63b4de8b4c9e08db6adc44bb1c8e1ef144 Mon Sep 17 00:00:00 2001 From: Olivier Sauter <Olivier.Sauter@epfl.ch> Date: Sat, 26 Oct 2019 07:48:04 +0200 Subject: [PATCH] add time_out as parameter to ask for specific time interval or instances for ids as well, at this stage dealt with before output of gdat_tcv for each case with single function gdat2time_out.m --- matlab/TCV/gdat_tcv.m | 340 ++++++++++++++++++-- matlab/TCV/tcv_help_parameters.m | 11 +- matlab/TCV/tcv_requests_mapping.m | 4 +- matlab/TCV_IMAS/tcv_get_ids_core_profiles.m | 117 ++++--- matlab/TCV_IMAS/tcv_get_ids_ec_antennas.m | 10 +- matlab/TCV_IMAS/tcv_get_ids_equilibrium.m | 197 ++++++++---- matlab/TCV_IMAS/tcv_get_ids_magnetics.m | 2 + matlab/TCV_IMAS/tcv_get_ids_nbi.m | 22 +- matlab/TCV_IMAS/tcv_get_ids_pf_active.m | 2 + matlab/TCV_IMAS/tcv_get_ids_tf.m | 5 +- matlab/TCV_IMAS/tcv_get_ids_wall.m | 4 +- matlab/TCV_IMAS/tcv_ids_bpol_probe.m | 14 +- matlab/TCV_IMAS/tcv_ids_circuit.m | 6 +- matlab/TCV_IMAS/tcv_ids_coil.m | 17 +- matlab/TCV_IMAS/tcv_ids_flux_loop.m | 7 +- matlab/TCV_IMAS/tcv_ids_headpart.m | 4 +- matlab/TCV_IMAS/tcv_ids_ip.m | 11 +- matlab/TCV_IMAS/tcv_ids_supply.m | 6 +- matlab/gdat2time_out.m | 288 +++++++++++++++++ 19 files changed, 891 insertions(+), 176 deletions(-) create mode 100644 matlab/gdat2time_out.m diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m index 9e781fde..fdeff219 100644 --- a/matlab/TCV/gdat_tcv.m +++ b/matlab/TCV/gdat_tcv.m @@ -181,7 +181,8 @@ if nargin>=2 && ivarargin_first_char~=1 return end data_request_eff = data_request.data_request; - data_request.data_request = lower(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 @@ -497,7 +498,14 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi') gdat_data.x = gdat_data.dim(dim_nontim); end end - if length(gdat_data.dim)>=mapping_for_tcv.timedim && mapping_for_tcv.timedim>0; gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim}; end + if length(gdat_data.dim)>=mapping_for_tcv.timedim && mapping_for_tcv.timedim>0; + gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim}; + mapping_for_tcv.gdat_timedim = mapping_for_tcv.timedim; + gdat_data.mapping_for.tcv = mapping_for_tcv; + end + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end gdat_data.units = aatmp.units; gdat_data.dimunits = aatmp.dimunits; if mapping_for_tcv.gdat_timedim>0 && mapping_for_tcv.gdat_timedim ~= mapping_for_tcv.timedim @@ -580,6 +588,9 @@ elseif strcmp(mapping_for_tcv.method,'expression') 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" @@ -687,6 +698,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.units = aatmp.units(end); gdat_data.dimunits = aatmp.dimunits; end + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'zgeom','z_geom'} @@ -721,6 +735,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if any(strcmp(fieldnames(zcontour),'units')) gdat_data.units = zcontour.units; end + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'b0'} @@ -766,6 +783,12 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.request_description = ['vacuum magnetic field at R0=' num2str(r0exp) 'm; COCOS=17']; gdat_data.r0 = r0exp; + % At this stage implement time_out at resulting structure as post processing, to allow various interpolation and other options. + % Should implement at server level at some point, could exend tdi_os + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'betan'} % 100*beta / |Ip[MA] * B0[T]| * a[m] @@ -788,6 +811,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.dim = beta.dim; gdat_data.t = beta.dim{1}; gdat_data.data = beta.data; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end ij=find(isfinite(ip.data)); ip_t = interp1(ip.dim{1}(ij),ip.data(ij),gdat_data.t); ij=find(isfinite(b0.data)); @@ -827,13 +853,13 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') source = [1 2 3]; end gdat_data.gdat_params.source = source; - if isfield(params_eff,'time_interval') && ~isempty(params_eff.time_interval) && length(params_eff.time_interval)>=2 - time_interval = params_eff.time_interval; + if isfield(params_eff,'cxrs_time_interval') && ~isempty(params_eff.cxrs_time_interval) && length(params_eff.cxrs_time_interval)>=2 + cxrs_time_interval = params_eff.cxrs_time_interval; cxrs_plot=1; else - time_interval = []; + cxrs_time_interval = []; end - gdat_data.gdat_params.time_interval = time_interval; + gdat_data.gdat_params.cxrs_time_interval = cxrs_time_interval; gdat_data.gdat_params.cxrs_plot = cxrs_plot; fit_tension_default = -100.; if isfield(params_eff,'fit_tension') @@ -859,7 +885,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') cxrs_params.prof.nc.taus = fit_tension.ni; cxrs_params.prof.zeff.taus = fit_tension.zeff; cxrs_params.k_plot = cxrs_plot; - cxrs_profiles = CXRS_get_profiles(shot,source,time_interval,cxrs_params); + cxrs_profiles = CXRS_get_profiles(shot,source,cxrs_time_interval,cxrs_params); inb_times = length(cxrs_profiles.Times); gdat_data.cxrs_params = cxrs_profiles.param; if isempty(cxrs_profiles.Times) || ~isfield(cxrs_profiles,'proffit') @@ -874,7 +900,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.(sub_eff_out).raw.error_bar = []; gdat_data.(sub_eff_out).raw.error_bar_rho = []; gdat_data.(sub_eff_out).raw.cxrs_system = []; - gdat_data.time_interval = []; + gdat_data.cxrs_time_interval = []; gdat_data.(sub_eff_out).units = sub_nodes_units{i}; if i==1 gdat_data.data = gdat_data.(sub_eff_out).fit.data; @@ -905,7 +931,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.(sub_eff_out).raw.error_bar = NaN * ones(inb_channels,inb_times); gdat_data.(sub_eff_out).raw.error_bar_rho = NaN * ones(inb_channels,inb_times); gdat_data.(sub_eff_out).raw.cxrs_system = NaN * ones(inb_channels,inb_times); - gdat_data.time_interval = []; + gdat_data.cxrs_time_interval = []; for it=1:inb_times if isfield(cxrs_profiles,sub_eff) nb_raw_points = length(cxrs_profiles.(sub_eff){it}.use.y); @@ -914,7 +940,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.(sub_eff_out).raw.error_bar(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.dy; gdat_data.(sub_eff_out).raw.error_bar_rho(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.dx; gdat_data.(sub_eff_out).raw.cxrs_system(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.sys; - gdat_data.time_interval{it} = cxrs_profiles.(sub_eff){it}.t_lim; + gdat_data.cxrs_time_interval{it} = cxrs_profiles.(sub_eff){it}.t_lim; end end gdat_data.(sub_eff_out).units = sub_nodes_units{i}; @@ -927,10 +953,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end gdat_data.t = cxrs_profiles.proffit.time; gdat_data.dim = {gdat_data.x; gdat_data.t}; - if isempty(time_interval) + if isempty(cxrs_time_interval) gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[' num2str(source) '],[],cxrs_params); % with cxrs_params']; else - gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[' num2str(source) '],[' num2str(time_interval) '],cxrs_params); % with cxrs_params']; + gdat_data.data_fullpath=['CXRS_get_profiles(' num2str(shot) ',[' num2str(source) '],[' num2str(cxrs_time_interval) '],cxrs_params); % with cxrs_params']; end gdat_data.dimunits{1} = ''; gdat_data.dimunits{2} = 's'; @@ -940,10 +966,48 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') else gdat_data = get_grids_1d(gdat_data,2,0,gdat_params.nverbose); end + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) && ~isempty(gdat_data.t) ... + && ~isempty(gdat_data.data) && isnumeric(gdat_data.data) + others_as_data={'error_bar'}; + size_data = size(gdat_data.data); + for i=1:length(sub_nodes_out) + fn_sub = fieldnames(gdat_data.(sub_nodes_out{i})); + fn_sub = setdiff(fn_sub,'units'); + for j=1:length(fn_sub) + fn_sub_sub=fieldnames(gdat_data.(sub_nodes_out{i}).(fn_sub{j})); + for k=1:length(fn_sub_sub) + if ~isempty(gdat_data.(sub_nodes_out{i}).(fn_sub{j}).(fn_sub_sub{k})) + try + size_leaf = size(gdat_data.(sub_nodes_out{i}).(fn_sub{j}).(fn_sub_sub{k})); + if size_data(gdat_data.mapping_for.tcv.timedim) == size_leaf(gdat_data.mapping_for.tcv.timedim) ... + && length(size_data)==length(size_leaf) + others_as_data{end+1} = [sub_nodes_out{i} '.' fn_sub{j} '.' fn_sub_sub{k}]; + end + catch + end + end + end + end + end + if gdat_data.gdat_params.nverbose >= 3 + disp(['in cxrs, fields with data reduced to match time_out:']) + others_as_data + end + [gdat_data] = gdat2time_out(gdat_data,1,others_as_data); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'eqdsk'} % + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time) + disp('time_out parameter not used for eqdsk, uses time') + gdat_data.gdat_params = rmfield(gdat_data.gdat_params,'time_out'); + else + warning('time_out parameter not relevant for eqdsk, use time to give time array for outputs') + return; + end + end time=1.; % default time if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time) time = gdat_data.gdat_params.time; @@ -1103,12 +1167,37 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') 'plot_eqdsk, write_eqdsk, read_eqdsk can be used']; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - case {'halpha'} - channels = -1; + case {'halphas'} + mapping.expression = '\base::pd:pd_001'; + channels = [1:18]; if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels) channels = gdat_data.gdat_params.channels; end % '\base::pd:pd_001'; + params_eff = gdat_data.gdat_params; + for i=1:length(channels) + params_eff.data_request=['\base::pd:pd_' num2str(channels(i),'%.3d')]; + try + halphas=gdat_tcv(gdat_data.shot,params_eff); % note: no need to set .doplot=0 since gdat_tcv does not call gdat_plot in any case + if i == 1; + gdat_data.data = NaN*ones(max(channels),length(halphas.t)); + gdat_data.t = halphas.t; + end + gdat_data.data(channels(i),:) = halphas.data; + catch + disp([params_eff.data_request ' not ok']); + end + end + gdat_data.label = ['\base::pd:pd_' num2str(min(channels),'%.3d') ' to ' num2str(max(channels),'%.3d')]; + gdat_data.x = channels; + gdat_data.dim{1} = gdat_data.x; + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits{1} = 'channel #'; + gdat_data.dimunits{2} = 's'; + gdat_data.mapping_for.tcv.gdat_timedim = 2; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'ids'} @@ -1251,6 +1340,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end end end + % create structure for icd sources from params and complete with defaults source_icd.ec = 'toray'; source_icd.nbi = ''; @@ -1263,7 +1353,6 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') source_icd.(gdat_data.gdat_params.source{i}) = gdat_data.gdat_params.(['source_' gdat_data.gdat_params.source{i}]); end end - mdsopen(shot); field_for_main_data = 'cd_tot'; % add each source in main.data, on ohm time array @@ -1280,6 +1369,13 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') else [pabs_gyro,icdtot,pow_dens,currentdrive_dens,rho_dep_pow,drho_pow,pmax,icdmax,currentdrive_dens_w2,rho_dep_icd,drho_icd] = astra_tcv_EC_exp(shot,[],[],[],[],[],gdat_data.gdat_params.trialindx); % centralized function for toray nodes end + if gdat_data.mapping_for.tcv.gdat_timedim ==2 + tgrid_to_change = {'pabs_gyro','icdtot','pow_dens','currentdrive_dens','rho_dep_pow','drho_pow','pmax', ... + 'icdmax','currentdrive_dens_w2','rho_dep_icd','drho_icd'}; + for i=1:length(tgrid_to_change) + eval([tgrid_to_change{i} '.tgrid = reshape(' tgrid_to_change{i} '.tgrid,1,numel(' tgrid_to_change{i} '.tgrid));']); + end + end ec_help = 'from toray icdint with extracting of effective Icd for given launcher depending on nb rays used'; % All EC related quantities, each substructure should have at least fields data,x,t,units,dim,dimunits,label to be copied onto gdat_data launchers_label = {'1','2','3','4','5','6','7','8','9','tot'}; @@ -1462,9 +1558,21 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ohm_help = ''; if strcmp(lower(source_icd.ohm),'ibs') ohm_data.cd_tot = gdat_tcv(shot,'\results::ibs:iohm'); + if gdat_data.mapping_for.tcv.gdat_timedim ==2 + tgrid_to_change = {'ohm_data.cd_tot.t','ohm_data.cd_tot.dim{1}'}; + for i=1:length(tgrid_to_change) + eval([tgrid_to_change{i} ' = reshape(' tgrid_to_change{i} ',1,numel(' tgrid_to_change{i} '));']); + end + end ohm_data.cd_tot.data = ohm_data.cd_tot.data'; ohm_data.cd_tot.units = strrep(ohm_data.cd_tot.units,'kA','A'); ohm_data.cd_dens = gdat_tcv(shot,'\results::ibs:johmav'); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R + if gdat_data.mapping_for.tcv.gdat_timedim ==2 + tgrid_to_change = {'ohm_data.cd_dens.t','ohm_data.cd_dens.dim{2}'}; + for i=1:length(tgrid_to_change) + eval([tgrid_to_change{i} ' = reshape(' tgrid_to_change{i} ',1,numel(' tgrid_to_change{i} '));']); + end + end ohm_data.cd_dens.units = strrep(ohm_data.cd_dens.units,'kA','A'); abc=get_grids_1d(ohm_data.cd_dens,1,1); ohm_data.cd_dens.rhotor_norm = abc.grids_1d.rhotornorm; @@ -1487,7 +1595,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if isempty(gdat_data.t); gdat_data.t = gdat_data.ohm.t; end gdat_data.ohm.data_fullpath = data_fullpath; gdat_data.ohm.help = ohm_help; - gdat_data.data(end+1,:) = interpos(gdat_data.ohm.t,gdat_data.ohm.data,gdat_data.t,-0.1); + ij_ok = [isfinite(gdat_data.ohm.data)]; + gdat_data.data(end+1,:) = interpos(-63,gdat_data.ohm.t(ij_ok),gdat_data.ohm.data(ij_ok),gdat_data.t,-0.1); gdat_data.label{end+1}=gdat_data.ohm.label; end % @@ -1496,9 +1605,21 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') bs_help = ''; if strcmp(lower(source_icd.bs),'ibs') bs_data.cd_tot = gdat_tcv(shot,'\results::ibs:ibs'); + if gdat_data.mapping_for.tcv.gdat_timedim ==2 + tgrid_to_change = {'bs_data.cd_tot.t','bs_data.cd_tot.dim{1}'}; + for i=1:length(tgrid_to_change) + eval([tgrid_to_change{i} ' = reshape(' tgrid_to_change{i} ',1,numel(' tgrid_to_change{i} '));']); + end + end bs_data.cd_tot.data = bs_data.cd_tot.data'; bs_data.cd_tot.units = strrep(bs_data.cd_tot.units,'kA','A'); bs_data.cd_dens = gdat_tcv(shot,'\results::ibs:jbsav'); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R + if gdat_data.mapping_for.tcv.gdat_timedim ==2 + tgrid_to_change = {'bs_data.cd_dens.t','bs_data.cd_dens.dim{2}'}; + for i=1:length(tgrid_to_change) + eval([tgrid_to_change{i} ' = reshape(' tgrid_to_change{i} ',1,numel(' tgrid_to_change{i} '));']); + end + end bs_data.cd_dens.units = strrep(bs_data.cd_dens.units,'kA','A'); abc=get_grids_1d(bs_data.cd_dens,1,1); bs_data.cd_dens.rhotor_norm = abc.grids_1d.rhotornorm; @@ -1521,7 +1642,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if isempty(gdat_data.t); gdat_data.t = gdat_data.bs.t; end gdat_data.bs.data_fullpath = data_fullpath; gdat_data.bs.help = bs_help; - gdat_data.data(end+1,:) = interpos(gdat_data.bs.t,gdat_data.bs.data,gdat_data.t,-0.1); + ij_ok = [isfinite(gdat_data.bs.data)]; + gdat_data.data(end+1,:) = interpos(-63,gdat_data.bs.t(ij_ok),gdat_data.bs.data(ij_ok),gdat_data.t,-0.1); gdat_data.label{end+1}=gdat_data.bs.label; end % @@ -1533,9 +1655,62 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.dim{2} = gdat_data.t; gdat_data.dimunits = {['index for each main source + total ' field_for_main_data], 's'}; gdat_data.data_fullpath = 'see in individual source substructure'; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + for i=1:length(gdat_data.gdat_params.source) + if ~isempty(gdat_data.(gdat_data.gdat_params.source{i}).data) && ~isempty(gdat_data.(gdat_data.gdat_params.source{i}).t) + % data with same time length as gdat_data.(gdat_data.gdat_params.source{i}).t + ab_tmp_rho = gdat_data.(gdat_data.gdat_params.source{i}); + ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for; + switch gdat_data.gdat_params.source{i} + case 'ec' + other_data = {'ec_data.p_abs_plasma.data','ec_data.p_abs_plasma.t', ... + ['ec_data.p_abs_plasma.dim{' num2str(gdat_data.mapping_for.tcv.gdat_timedim) '}'], ... + 'ec_data.cd_tot.data','ec_data.cd_tot.t', ... + ['ec_data.cd_tot.dim{' num2str(gdat_data.mapping_for.tcv.gdat_timedim) '}']}; + case 'ohm' + other_data = {'ohm_data.cd_tot.data','ohm_data.cd_tot.t', ['ohm_data.cd_tot.dim{1}']}; + case 'bs' + other_data = {'bs_data.cd_tot.data','bs_data.cd_tot.t', ['bs_data.cd_tot.dim{1}']}; + otherwise + disp(gdat_data.gdat_params.source{i}) + keyboard + end + [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,other_data); + gdat_data.(gdat_data.gdat_params.source{i}) = rmfield(ab_tmp_rho,{'gdat_params','mapping_for'}); + end + end + % extra ec subfields + if ~isempty(gdat_data.ec.data) && ~isempty(gdat_data.ec.t) && isfield(gdat_data.ec,'ec_data') + extra_fields=fieldnames(gdat_data.ec.ec_data); + extra_fields = setdiff(extra_fields,{'p_abs_plasma','cd_tot'}); + for i=1:length(extra_fields) + ab_tmp_rho = gdat_data.ec.ec_data.(extra_fields{i}); + ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for; + ab_tmp_rho.mapping_for.tcv.gdat_timedim = length(size(ab_tmp_rho.data)); % time always last dim for these fields + [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1); + gdat_data.ec.ec_data.(extra_fields{i}) = rmfield(ab_tmp_rho,{'gdat_params','mapping_for'}); + end + end + % extra cd_dens fields in ohm, bs + extra_fields = intersect({'bs','ohm'},gdat_data.gdat_params.source); + for i=1:length(extra_fields) + if ~isempty(gdat_data.(extra_fields{i}).data) && ~isempty(gdat_data.(extra_fields{i}).t) && ... + isfield(gdat_data.(extra_fields{i}),[extra_fields{i} '_data']) + gdat_data.(extra_fields{i}).([extra_fields{i} '_data']).cd_dens.gdat_params.time_out = gdat_data.gdat_params.time_out; + gdat_data.(extra_fields{i}).([extra_fields{i} '_data']).cd_dens = gdat2time_out( ... + gdat_data.(extra_fields{i}).([extra_fields{i} '_data']).cd_dens,1,{'rhotor_norm'}); + end + end + + [gdat_data] = gdat2time_out(gdat_data,1); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'mhd'} + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + warning(['parameter time_out not implemented yet, will load full time data; ask O. Sauter if needed' char(10)]); + gdat_data.gdat_params = rmfield(gdat_data.gdat_params,'time_out'); + end if isfield(gdat_data.gdat_params,'nfft') && ~isempty(gdat_data.gdat_params.nfft) % used in gdat_plot for spectrogram plot else @@ -1628,6 +1803,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.gdat_params.edge = 0; end [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.edge,gdat_params.nverbose); + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1,{'error_bar'}); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'ne_rho', 'te_rho', 'nete_rho'} @@ -1662,7 +1840,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') else % calculate new psiscatvol [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, edge_str_dot, psitbx_str, gdat_params); - + % Add the results to the output of gdat gdat_data.psiscatvol = psiscatvol; gdat_data.psi_max = psi_max; @@ -1686,7 +1864,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.dim{1}=rho; gdat_data.dimunits=[{'sqrt(psi_norm)'} ; {'time [s]'}]; gdat_data.x=rho; - % add grids_1d to have rhotor, etc if cxrs_rho was asked for + % add grids_1d to have rhotor, etc gdat_data = get_grids_1d(gdat_data,2,1,gdat_params.nverbose); %%%%%%%%%%% @@ -1876,6 +2054,40 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') else gdat_data.gdat_params.fit_type = default_fit_type; end + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + % add extra fields then correct + ab_tmp_rho = gdat_data.fit; ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for; + [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,{'rhotornorm','rhovolnorm'}); + gdat_data.fit = rmfield(ab_tmp_rho,{'gdat_params','mapping_for'}); + % + gdat_data.psi_max.dim{2} = gdat_data.psi_max.dim{1}; + extra_fields = {'error_bar','x','dim{1}','psiscatvol.data','psiscatvol.dim','psi_max.data','psi_max.dim'}; + if strcmp(data_request_eff,'ne_rho') + extra_fields(end+1:end+3) = {'firrat','data_raw','error_bar_raw'}; + gdat_data.firrat = reshape(gdat_data.firrat,1,numel(gdat_data.firrat)); + end + [gdat_data] = gdat2time_out(gdat_data,22,extra_fields); + gdat_data.psi_max.dim = gdat_data.psi_max.dim{2}; + if strcmp(data_request_eff,'nete_rho') + extra_sub = {'ne','te'}; + for i=1:length(extra_sub) + ab_tmp_rho = gdat_data.fit.(extra_sub{i}); ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for; + [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,{'rhotornorm','rhovolnorm'}); + gdat_data.fit.(extra_sub{i}) = rmfield(ab_tmp_rho,{'gdat_params','mapping_for'}); + % + ab_tmp_rho = gdat_data.(extra_sub{i}); + ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for; + ab_tmp_rho.t = gdat_data.t; + extra_fields = {'error_bar'}; + if strcmp(extra_sub{i},'ne'); + extra_fields(end+1:end+3) = {'firrat','data_raw','error_bar_raw'}; + ab_tmp_rho.firrat = reshape(ab_tmp_rho.firrat,1,numel(ab_tmp_rho.firrat)); + end + [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,extra_fields); + gdat_data.(extra_sub{i}) = rmfield(ab_tmp_rho,{'gdat_params','mapping_for','t'}); + end + end + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'powers'} @@ -1944,7 +2156,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.dimunits{1} = 's'; gdat_data.units = 'W'; tension = -1e5; - vloop_smooth=interpos(vloop.t,vloop.data,gdat_data.ohm.t,tension); + vloop_smooth=interpos(-63,vloop.t,vloop.data,gdat_data.ohm.t,tension); ip_t = interp1(ip.t,ip.data,gdat_data.ohm.t); gdat_data.ohm.data = -vloop_smooth.*ip_t; % TCV has wrong sign for Vloop gdat_data.ohm.raw_data = -vloop.data.*ip_t; @@ -2063,6 +2275,19 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.label{end+1}='total heating'; gdat_data.dim{2} = gdat_data.x; gdat_data.dimunits = {'s', 'index for each source + total heating'}; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + for i=1:length(sources_avail) + ab_tmp_rho = gdat_data.(sources_avail{i}); + ab_tmp_rho.gdat_params = gdat_data.gdat_params; ab_tmp_rho.mapping_for = gdat_data.mapping_for; + if strcmp(sources_avail{i},'nbi') + [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,{'energy'}); + else + [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1); + end + gdat_data.(sources_avail{i}) = rmfield(ab_tmp_rho,{'gdat_params','mapping_for'}); + end + [gdat_data] = gdat2time_out(gdat_data,1); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'q_rho'} @@ -2097,6 +2322,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.request_description = 'q(rhopol\_norm)'; % add grids_1d to have rhotor, etc gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose); + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,21); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'rbphi_rho', 'rbtor_rho'} @@ -2133,6 +2361,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.request_description = nodenameeff; % add grids_1d to have rhotor, etc gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose); + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + % add extra fields then correct + [gdat_data] = gdat2time_out(gdat_data,21); + end case {'phi_tor', 'phitor', 'toroidal_flux'} % Phi(LCFS) = int(Bphi dSphi), can use Eq.(11) of "Tokamak coordinate conventions: COCOS" paper: @@ -2231,6 +2463,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.units = phi_tor.units; if (length(gdat_data.dim)>=mapping_for_tcv.gdat_timedim); gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim}; end if (length(gdat_data.dim)>0); gdat_data.x = gdat_data.dim{1}; end + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + % add extra fields then correct + [gdat_data] = gdat2time_out(gdat_data,1); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'pprime', 'pprime_rho'} @@ -2255,8 +2491,11 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.dimunits = tracetdi.dimunits; gdat_data.units = tracetdi.units; gdat_data.request_description = 'pprime=dp/dpsi'; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end - case {'pressure', 'pressure_rho'} + case {'pressure', 'pressure_rho', 'p_rho'} if liuqe_matlab==0 nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('p_rho',1,0) '",''' psitbx_str ''')']; tracetdi=tdi(nodenameeff); @@ -2278,11 +2517,15 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.dimunits = tracetdi.dimunits; gdat_data.units = tracetdi.units; gdat_data.request_description = 'pressure'; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'psi_edge'} % psi at edge, 0 by construction in Liuqe, thus not given % add surface_psi from surface_flux and d(surface_flux)/dt = vloop +keyboard nodenameeff=['\results::psi_axis' substr_liuqe]; if liuqe_version_eff==-1 nodenameeff=[begstr 'q_psi' substr_liuqe]; @@ -2308,6 +2551,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.surface_psi = surface_psi.data; gdat_data.vsurf = vsurf; gdat_data.vsurf_description = ['time derivative from surface_psi, obtained from \results::surface_flux']; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'r_contour_edge', 'z_contour_edge'} @@ -2327,6 +2573,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.data_fullpath=nodenameeff; gdat_data.dimunits = tracetdi.dimunits; gdat_data.units = tracetdi.units; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'rhotor_edge', 'rhotor', 'rhotor_norm'} @@ -2426,7 +2675,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if strcmp(data_request_eff,'rhotor_edge') gdat_data.data = gdat_data.rhotor_edge; gdat_data.dim = phi_tor.dim(2); % while there is a problem with \tcv_shot::top.results.equil... dim nodes - gdat_data.dimunits = phi_tor.dimunits{2}; + gdat_data.dimunits = phi_tor.dimunits(2); gdat_data.units = 'm'; gdat_data.data_fullpath=['rhotor_edge=sqrt(Phi(edge)/pi/B0) from phi_tor']; elseif strcmp(data_request_eff,'rhotor_norm') @@ -2452,6 +2701,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.rhotor_edge = []; end end + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + % nothing to do since above uses gdat call getting reduced time already + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'rhovol','rho_vol','volume_rho','volume'} @@ -2484,15 +2736,15 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end gdat_data.units = tracetdi.units; if strcmp(data_request_eff,'volume') - if length(gdat_data.dim >=2) - gdat_data.data = tracetdi.data(end,:); + if length(tracetdi.dim) >= 2 + gdat_data.data = reshape(tracetdi.data(end,:),numel(tracetdi.data(end,:)),1); % needed consistent with timedim for gdat2time gdat_data.volume_edge = gdat_data.data; - gdat_data.dim{1} = tracetdi.dim{2}; + gdat_data.dim{1} = reshape(tracetdi.dim{2},numel(tracetdi.dim{2}),1); gdat_data.dimunits{1} = tracetdi.dimunits{2}; else mapping_for_tcv.gdat_timedim = 1; - gdat_data.data = tracetdi.data; - gdat_data.dim = tracetdi.dim; + gdat_data.data = reshape(tracetdi.data,numel(tracetdi.data),1); + gdat_data.dim = reshape(tracetdi.dim,numel(tracetdi.dim),1); gdat_data.dimunits = tracetdi.dimunits; gdat_data.volume_edge = gdat_data.data; end @@ -2520,6 +2772,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end end gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim}; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1,{'volume_edge'}); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'rtc'} @@ -2714,8 +2969,20 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end return end - if ~isfield(gdat_data.gdat_params,'time_interval') - gdat_data.gdat_params.time_interval = []; + time_interval = []; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + if length(gdat_data.gdat_params.time_out) == 2 + time_interval = gdat_data.gdat_params.time_out; + else + if length(gdat_data.gdat_params.time_out) == 1 + % 200ms includes all characteristic time constants + time_interval = [gdat_data.gdat_params.time_out-0.1 gdat_data.gdat_params.time_out+0.1]; + else + time_interval = [min(gdat_data.gdat_params.time_out)-0.1 max(gdat_data.gdat_params.time_out)+0.1]; + end + warning(['Expects a time interval [t1 t2] for ' data_request_eff ' in param time_out, uses [' ... + num2str(time_interval(1)) ',' num2str(time_interval(2)) ']' char(10)]) + end end if ~isfield(gdat_data.gdat_params,'freq') gdat_data.gdat_params.freq = 'slow'; @@ -2737,7 +3004,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') freq_opt = 0; if strcmp(gdat_data.gdat_params.freq,'fast'); freq_opt = 1; end t_int = [0 10]; % starts from 0 otherwise mpxdata gives data from t<0 - if ~isempty(gdat_data.gdat_params.time_interval); t_int = gdat_data.gdat_params.time_interval; end + if ~isempty(time_interval); t_int = time_interval; end gdat_data.top.data = []; gdat_data.top.x = []; gdat_data.top.channel = []; @@ -2814,11 +3081,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end try - if isempty(gdat_data.gdat_params.time_interval); + if isempty(time_interval); [sig,t,Sat_channel,cat_default] = XTOMOGetData(shot,0,4,camera_xtomo,[]); else - [sig,t,Sat_channel,cat_default] = XTOMOGetData(shot,gdat_data.gdat_params.time_interval(1), ... - gdat_data.gdat_params.time_interval(2),camera_xtomo,[]); + [sig,t,Sat_channel,cat_default] = XTOMOGetData(shot,time_interval(1),time_interval(2),camera_xtomo,[]); end catch if gdat_data.gdat_params.nverbose>=1 @@ -2877,6 +3143,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.dimunits = tracetdi.dimunits; gdat_data.units = tracetdi.units; gdat_data.request_description = 'ttprime=t*dt/dpsi'; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'profnerho','profterho'} @@ -2929,6 +3198,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end gdat_data.request_description = 'quick autofits within thomson nodes, using version'; gdat_data.fullpath = ['Thomson autfits from ' nodenameeff]; + gdat_data.mapping_for.tcv.gdat_timedim = 2; + if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out) + [gdat_data] = gdat2time_out(gdat_data,1); + end otherwise if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_tcv']); end @@ -3169,4 +3442,3 @@ psiscatvol.dim{2} = t_th; psi_norm = psitbxp2p(psi,'01'); psi_max.data = psi_norm.psimag(i_psi); psi_max.dim = {t_th}; - diff --git a/matlab/TCV/tcv_help_parameters.m b/matlab/TCV/tcv_help_parameters.m index 670d2293..289e15be 100644 --- a/matlab/TCV/tcv_help_parameters.m +++ b/matlab/TCV/tcv_help_parameters.m @@ -5,9 +5,9 @@ function help_struct = tcv_help_parameters(parameter_list) % % return the whole help structure if parameter_list empty or not provided % -% do: +% do: % help_struct = tcv_help_parameters(fieldnames(gdat_data.gdat_params)); -% +% % to get relevant help description % @@ -24,13 +24,14 @@ help_struct_all = struct(... % TCV related help_struct_all.cxrs_plot = '0 (default) no plots, 1 get plots from CXRS_get_profiles see ''help CXRS_get_profiles'' for k_plot values'; -help_struct_all.time_interval = ['if provided sets a specific time interval [tstart tend].' ... - char(10) 'cxrs: (time_interval can have several nbs) take data and average over time interval(s) only, plots from CXRS_get_profiles are then provided' ... +help_struct_all.cxrs_time_interval = ['cxrs: (time_interval can have several nbs) take data and average over time interval(s) only, plots from CXRS_get_profiles are then provided' ... ' as well']; help_struct_all.fit_tension = ['smoothing value used in interpos fitting routine, -30 means ''30 times default value'', thus -1 often a' ... ' good value' char(10) ... 'cxrs: if numeric, default for all cases, if structure, default for non given fields']; -help_struct_all.time = 'time(s) value(s) if relevant, for example eqdsk is provided by default only for time=1.0s'; +help_struct_all.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.zshift = 'vertical shift of equilibrium, either for eqdsk (1 to shift to zaxis=0) or for mapping measurements on to rho surfaces [m]'; help_struct_all.cocos = ['cocos value desired in output, uses eqdsk_cocos_transform. Note should use latter if a specific Ip and/or B0 sign' ... 'is wanted. See O. Sauter et al Comput. Phys. Commun. 184 (2013) 293']; diff --git a/matlab/TCV/tcv_requests_mapping.m b/matlab/TCV/tcv_requests_mapping.m index 3d64162f..b503f245 100644 --- a/matlab/TCV/tcv_requests_mapping.m +++ b/matlab/TCV/tcv_requests_mapping.m @@ -245,7 +245,7 @@ switch lower(data_request) 'tmp_data2=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ... 'gdat_tmp.data = gdat_tmp.data./(tmp_data2+1e-5);']; case {'ec_data', 'aux', 'h_cd', 'nbi_data', 'ic_data', 'lh_data', 'ohm_data', 'bs_data'} - mapping.timedim = 1; + mapping.timedim = 2; mapping.label = 'various Pdens, Icd, jcd'; mapping.method = 'switchcase'; case {'phi_tor', 'phitor', 'toroidal_flux'} @@ -288,7 +288,7 @@ switch lower(data_request) mapping.timedim = 2; mapping.label = 'pprime'; mapping.method = 'switchcase'; - case {'pressure', 'p_rho'} % note: not pressure from liuqe fortran which is 2D + case {'pressure', 'p_rho', 'pressure_rho'} % note: not pressure from liuqe fortran which is 2D mapping.timedim = 2; mapping.label = 'pressure'; mapping.method = 'switchcase'; diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m b/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m index 02a68e75..5d497abe 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m +++ b/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m @@ -25,12 +25,22 @@ if exist('gdat_params') [ids_core_profiles, params_cores_profiles] = tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles','gdat_params',gdat_params,varargin{:}); else [ids_core_profiles, params_cores_profiles] = tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles',varargin{:}); + aa=gdat_tcv; + gdat_params = aa.gdat_params; % to get default params end +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) + ids_core_profiles_description = []; % base all from times fro nete_rho (which should be conf by default) -temp_1d.ne_rho = gdat(params_cores_profiles.shot,'ne_rho','machine',machine,'fit',1); -temp_1d.te_rho = gdat(params_cores_profiles.shot,'te_rho','machine',machine,'fit',1); +params_eff = params_eff_ref; +params_eff.data_request='ne_rho'; params_eff.fit = 1; +temp_1d.ne_rho = gdat(params_cores_profiles.shot,params_eff); +temp_1d_desc.ne_rho = params_eff.data_request; +params_eff.data_request = 'te_rho'; +temp_1d.te_rho = gdat(params_cores_profiles.shot,params_eff); +temp_1d_desc.te_rho = params_eff.data_request; % compute grids_1d for fit rhopol grid, need a gdat_data like structure temp_1d.fit.te_rho = temp_1d.te_rho.fit; temp_1d.fit.te_rho.gdat_params = temp_1d.te_rho.gdat_params; @@ -60,55 +70,67 @@ ids_core_profiles.profiles_1d(1:length(ids_core_profiles.time)) = ids_core_profi % % vacuum_toroidal_field and time, using homogeneous % -vacuum_toroidal_field.b0=gdat(params_cores_profiles.shot,'b0','source','liuqe','machine',machine); % to get on liuqe time array -vacuum_toroidal_field_desc.b0 = '''b0'',''source'',''liuqe'''; -vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'',''source'',''liuqe'']'; +params_eff = params_eff_ref; +params_eff.data_request='b0'; +vacuum_toroidal_field.b0=gdat(params_cores_profiles.shot,params_eff); +vacuum_toroidal_field_desc.b0 = params_eff.data_request; ids_core_profiles.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0; ids_core_profiles.vacuum_toroidal_field.b0 = interpos(63,vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data,ids_core_profiles.time,-1); -ids_core_profiles_description.vacuum_toroidal_field = vacuum_toroidal_field_desc; +ids_core_profiles_description.vacuum_toroidal_field = [vacuum_toroidal_field_desc.b0 ' on ids_core_profiles.time, with interpos(63)']; % global_quantities data into local global_quantities.* structure with correct end names and global_quantities_desc.* with description. Use temp.* and temp_desc.* structures for temporary data +params_eff.data_request='ip'; global_quantities.ip = gdat(params_cores_profiles.shot,'ip','machine',machine); -global_quantities_desc.ip = 'ip'; +global_quantities_desc.ip = params_eff.data_request; -current_non_inductive = gdat(params_cores_profiles.shot,'ohm_data','machine',machine); +params_eff.data_request = 'ohm_data'; +current_non_inductive = gdat(params_cores_profiles.shot,params_eff); global_quantities.current_non_inductive.data = current_non_inductive.ohm.ohm_data.cd_tot.data; global_quantities.current_non_inductive.t = current_non_inductive.ohm.ohm_data.cd_tot.t; global_quantities_desc.current_non_inductive = current_non_inductive.ohm.ohm_data.cd_tot.label; -current_bootstrap = gdat(params_cores_profiles.shot,'bs_data','machine',machine); +params_eff.data_request = 'bs_data'; +current_bootstrap = gdat(params_cores_profiles.shot,params_eff); global_quantities.current_bootstrap.data = current_bootstrap.bs.bs_data.cd_tot.data; global_quantities.current_bootstrap.t = current_bootstrap.bs.bs_data.cd_tot.t; global_quantities_desc.current_bootstrap = current_bootstrap.bs.bs_data.cd_tot.label; -global_quantities.v_loop = gdat(params_cores_profiles.shot,'vloop','machine',machine); -global_quantities_desc.v_loop = 'vloop'; -global_quantities.li_3 = gdat(params_cores_profiles.shot,'li','machine',machine); -global_quantities_desc.li_3 = 'li'; -global_quantities.beta_tor = gdat(params_cores_profiles.shot,'beta','machine',machine); -global_quantities_desc.beta_tor = 'beta'; -global_quantities.beta_tor_norm = gdat(params_cores_profiles.shot,'betan','machine',machine); -global_quantities_desc.beta_tor_norm = 'betan'; -global_quantities.beta_pol = gdat(params_cores_profiles.shot,'betap','machine',machine); -global_quantities_desc.beta_pol = 'betap'; -global_quantities.energy_diamagnetic = gdat(params_cores_profiles.shot,'w_mhd','machine',machine); -global_quantities_desc.energy_diamagnetic = 'w_mhd'; -global_quantities.z_eff_resistive = gdat(params_cores_profiles.shot,'results.conf:z_eff','machine',machine,'fit',1); -global_quantities_desc.z_eff_resistive = 'results.conf:z_eff'; +params_eff.data_request = 'vloop'; +global_quantities.v_loop = gdat(params_cores_profiles.shot,params_eff); +global_quantities_desc.v_loop = params_eff.data_request; +params_eff.data_request = 'li'; +global_quantities.li_3 = gdat(params_cores_profiles.shot,params_eff); +global_quantities_desc.li_3 = params_eff.data_request; +params_eff.data_request = 'beta'; +global_quantities.beta_tor = gdat(params_cores_profiles.shot,params_eff); +global_quantities_desc.beta_tor = params_eff.data_request; +params_eff.data_request = 'betan'; +global_quantities.beta_tor_norm = gdat(params_cores_profiles.shot,params_eff); +global_quantities_desc.beta_tor_norm = params_eff.data_request; +params_eff.data_request = 'betap'; +global_quantities.beta_pol = gdat(params_cores_profiles.shot,params_eff); +global_quantities_desc.beta_pol = params_eff.data_request; +params_eff.data_request = 'w_mhd'; +global_quantities.energy_diamagnetic = gdat(params_cores_profiles.shot,params_eff); +global_quantities_desc.energy_diamagnetic = params_eff.data_request; +params_eff.data_request = 'results.conf:z_eff'; params_eff_fit1=params_eff; params_eff_fit1.fit=1; +global_quantities.z_eff_resistive = gdat(params_cores_profiles.shot,params_eff_fit1); +global_quantities_desc.z_eff_resistive = params_eff_fit1.data_request; ids_core_profiles_description.global_quantities = global_quantities_desc; - global_quantities_fieldnames = fieldnames(global_quantities); +ids_core_profiles_global_quantities_fieldnames = fieldnames(ids_core_profiles.global_quantities); +global_quantities_fieldnames_eff = intersect(global_quantities_fieldnames,ids_core_profiles_global_quantities_fieldnames); special_fields = {''}; % fields needing non-automatic treatments -for i=1:length(global_quantities_fieldnames) - if ~any(strcmp(global_quantities_fieldnames{i},special_fields)) - if ~isstruct(ids_core_profiles.global_quantities.(global_quantities_fieldnames{i})) - ids_core_profiles.global_quantities.(global_quantities_fieldnames{i}) = ... - interpos(global_quantities.(global_quantities_fieldnames{i}).t, ... - global_quantities.(global_quantities_fieldnames{i}).data,ids_core_profiles.time,tens_time); +for i=1:length(global_quantities_fieldnames_eff) + if ~any(strcmp(global_quantities_fieldnames_eff{i},special_fields)) + if ~isstruct(ids_core_profiles.global_quantities.(global_quantities_fieldnames_eff{i})) + ids_core_profiles.global_quantities.(global_quantities_fieldnames_eff{i}) = ... + interpos(global_quantities.(global_quantities_fieldnames_eff{i}).t, ... + global_quantities.(global_quantities_fieldnames_eff{i}).data,ids_core_profiles.time,tens_time); else - special_fields{end+1} = global_quantities_fieldnames{i}; + special_fields{end+1} = global_quantities_fieldnames_eff{i}; end end end @@ -117,11 +139,15 @@ end % %% profiles_1d (cannot use eqdsk since not same radial mesh) -temp_1d.area = gdat(shot,'area_rho','machine',machine); +params_eff.data_request = 'area_rho'; +temp_1d.area = gdat(params_cores_profiles.shot,params_eff); +temp_1d_desc.area = params_eff.data_request; for ir=1:length(temp_1d.area.x) area_cpt(ir,:) = interpos(temp_1d.area.t,temp_1d.area.data(ir,:),ids_core_profiles.time,tens_time); end -temp_1d.q = gdat(shot,'q_rho','machine',machine); +params_eff.data_request = 'q_rho'; +temp_1d.q = gdat(params_cores_profiles.shot,params_eff); +temp_1d_desc.q = params_eff.data_request; for ir=1:length(temp_1d.q.x) q_cpt(ir,:) = interpos(temp_1d.q.t,temp_1d.q.data(ir,:),ids_core_profiles.time,tens_time); end @@ -196,9 +222,15 @@ switch error_bar end % ion, assume only D if no CXRS (need to ask how to check if H...) -temp_1d.cxrs_rho = gdat(params_cores_profiles.shot,'cxrs','machine',machine,'fit',1); -temp_1d.ti_conf_rho = gdat(params_cores_profiles.shot,'results.conf:ti','machine',machine,'fit',1); -temp_1d.ni_conf_rho = gdat(params_cores_profiles.shot,'results.conf:ni','machine',machine,'fit',1); +params_eff_fit1.data_request = 'cxrs'; +temp_1d.cxrs_rho = gdat(params_cores_profiles.shot,params_eff_fit1); +temp_1d_desc.cxrs_rho = params_eff_fit1.data_request; +params_eff_fit1.data_request = 'results.conf:ti'; +temp_1d.ti_conf_rho = gdat(params_cores_profiles.shot,params_eff_fit1); +temp_1d_desc.ti_conf_rho = params_eff_fit1.data_request; +params_eff_fit1.data_request = 'results.conf:ni'; +temp_1d.ni_conf_rho = gdat(params_cores_profiles.shot,params_eff_fit1); +temp_1d_desc.ni_conf_rho = params_eff_fit1.data_request; data_fullpath_raw = 'Ti(C sometimes B) from cxrs system 1 to 3'; temp_1d.ti.raw = temp_1d.cxrs_rho.ti.raw; temp_1d.ti.raw.shot = temp_1d.cxrs_rho.shot;temp_1d.ti.raw.gdat_params = temp_1d.cxrs_rho.gdat_params; @@ -209,10 +241,13 @@ if ~isempty(temp_1d.cxrs_rho.ti.raw.data) else data_fullpath_fit = 'Ti from fit in CONF node, using Te/Ti and equilibrium Wmhd'; end +temp_1d_desc.ti.raw = data_fullpath_fit; temp_1d.ti.fit = temp_1d.ti_conf_rho; temp_1d.ti.fit =get_grids_1d(temp_1d.ti.fit,1,1); +temp_1d_desc.ti.fit = temp_1d_desc.ti_conf_rho; temp_1d.ni.fit = temp_1d.ni_conf_rho; temp_1d.ni.fit =get_grids_1d(temp_1d.ni.fit,1,1); +temp_1d_desc.ni.fit = temp_1d_desc.ni_conf_rho; it_ti = iround_os(temp_1d.ti.fit.t,ids_core_profiles.time); % assumed 1 impurity with Zp=6 Zp = 6.; @@ -268,7 +303,6 @@ if ~isempty(temp_1d.cxrs_rho.ti.fit.data) end end -temp_1d.q = gdat(shot,'q_rho','machine',machine); for ir=1:length(temp_1d.q.x) q_cpt(ir,:) = interpos(temp_1d.q.t,temp_1d.q.data(ir,:),ids_core_profiles.time,tens_time); end @@ -289,11 +323,16 @@ for it=1:length(ids_core_profiles.time) end end +ids_core_profiles_description.temp_1d_desc = temp_1d_desc; +ids_core_profiles_description.profiles_1d.magnetic_shear = 'from interpos with rhotor'; + if nargin <= 2 - ids_core_profiles.code.name = ['tcv_get_ids_core_profiles, within gdat, with shot= ' num2str(shot) ]; + ids_core_profiles.code.name = ['tcv_get_ids_core_profiles, within gdat, with shot= ' num2str(params_cores_profiles.shot) ]; else - ids_core_profiles.code.name = ['tcv_get_ids_core_profiles, within gdat, with shot= ' num2str(shot) '; varargin: ' varargin{:}]; + ids_core_profiles.code.name = ['tcv_get_ids_core_profiles, within gdat, with shot= ' num2str(params_cores_profiles.shot) '; varargin: ' varargin{:}]; end +ids_core_profiles_description.code.name = ids_core_profiles.code.name; + ids_core_profiles.code.output_flag = zeros(size(ids_core_profiles.time)); % make arrays not filled in empty: diff --git a/matlab/TCV_IMAS/tcv_get_ids_ec_antennas.m b/matlab/TCV_IMAS/tcv_get_ids_ec_antennas.m index 717c393d..4985b651 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_ec_antennas.m +++ b/matlab/TCV_IMAS/tcv_get_ids_ec_antennas.m @@ -11,15 +11,21 @@ if exist('gdat_params') 'gdat_params',gdat_params,varargin{:}); else [ids_ec_antennas, params_ec_antennas] = tcv_ids_headpart(shot,ids_ec_antennas_empty,'ec_antennas','homogeneous_time',0,varargin{:}); + aa=gdat_tcv; + gdat_params = aa.gdat_params; % to get default params end +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) % As a general rule, for a new substructure under the main ids, construct a local structure like: % "global_quantities" with subfields being the relevant data to get and a local structure: % "global_quantities_desc" which contains the same subfields themselves containing the gdat string aftre shot used % -pow=gdat_tcv(shot,'powers'); - +params_eff = params_eff_ref; +params_eff.data_request = 'powers'; +pow=gdat_tcv(shot,params_eff); +pow_desc = params_eff.data_request; if isempty(pow.ec.t) || isempty(pow.ec.data) || max(pow.ec.data(:,end))<1e-10 ids_ec_antennas_description.comment = 'no power'; return diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m index 23a5ef28..195391d1 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m +++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m @@ -10,6 +10,8 @@ if exist('gdat_params') [ids_equilibrium, params_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium','gdat_params',gdat_params,varargin{:}); else [ids_equilibrium, params_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium',varargin{:}); + aa=gdat_tcv; + gdat_params = aa.gdat_params; % to get default params end % As a general rule, for a new substructure under the main ids, construct a local structure like: @@ -19,8 +21,14 @@ end % vacuum_toroidal_field and time, using homogeneous % %% liuqe.m at this stage is missing correction 0.996, so use std source by time of default liuqe to make sure -bb = gdat(params_equilibrium.shot,'b0','source','liuqe','machine',gdat_params.machine,'liuqe',gdat_params.liuqe); % to get liuqe time array -vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,'b0','machine',gdat_params.machine); +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) + +params_eff = params_eff_ref; +params_eff.data_request='b0'; params_eff.source='liuqe'; % to get liuqe time array +bb = gdat(params_equilibrium.shot,params_eff); +params_eff = rmfield(params_eff,'source'); % to get original magnetics data +vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,params_eff); vacuum_toroidal_field.b0.data = interp1(vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data,bb.t); vacuum_toroidal_field.b0.t = bb.t; vacuum_toroidal_field.b0.dim = {vacuum_toroidal_field.b0.t}; @@ -46,40 +54,57 @@ ids_equilibrium.time_slice(1:numel(ids_equilibrium.time)) = ids_equilibrium.time % $$$ end % $$$ temp_desc.eqdsks{1} = '''eqdsk'',''time'',ids_equilibrium.time(it)'; -global_quantities.area = gdat(params_equilibrium.shot,'area_edge','machine',gdat_params.machine); -global_quantities_desc.area = 'area_edge'; -global_quantities.beta_normal = gdat(params_equilibrium.shot,'betan','machine',gdat_params.machine); -global_quantities_desc.beta_normal = 'betan'; -global_quantities.beta_pol = gdat(params_equilibrium.shot,'betap','machine',gdat_params.machine); -global_quantities_desc.beta_pol = 'betap'; -global_quantities.beta_tor = gdat(params_equilibrium.shot,'beta','machine',gdat_params.machine); -global_quantities_desc.beta_tor = 'beta'; -global_quantities.energy_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',gdat_params.machine); -global_quantities_desc.energy_mhd = 'w_mhd'; -global_quantities.ip = gdat(params_equilibrium.shot,'ip','machine',gdat_params.machine); -global_quantities_desc.ip = 'ip'; +params_eff = params_eff_ref; +params_eff.data_request = 'area_edge'; +global_quantities.area = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.area = params_eff.data_request; +params_eff.data_request = 'betan'; +global_quantities.beta_normal = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.beta_normal = params_eff.data_request; +params_eff.data_request = 'betap'; +global_quantities.beta_pol = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.beta_pol = params_eff.data_request; +params_eff.data_request = 'beta'; +global_quantities.beta_tor = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.beta_tor = params_eff.data_request; +params_eff.data_request = 'w_mhd'; +global_quantities.energy_mhd = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.energy_mhd = params_eff.data_request; +params_eff.data_request = 'ip'; +global_quantities.ip = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.ip = params_eff.data_request; % length_pol = gdat(params_equilibrium.shot,'length_pol','machine',gdat_params.machine); % to be added -global_quantities.li_3 = gdat(params_equilibrium.shot,'li','machine',gdat_params.machine); -global_quantities_desc.li_3 = 'li'; -temp.r_magnetic_axis = gdat(params_equilibrium.shot,'r_axis','machine',gdat_params.machine); -temp_desc.r_magnetic_axis = 'r_axis'; -temp.z_magnetic_axis = gdat(params_equilibrium.shot,'z_axis','machine',gdat_params.machine); -temp_desc.z_magnetic_axis = 'z_axis'; -temp.psi_axis = gdat(params_equilibrium.shot,'psi_axis','machine',gdat_params.machine); % needs to add psi_edge sincepsi_axis liuqe assuming 0 dege value -temp_desc.psi_axis = 'psi_axis'; -global_quantities.psi_boundary = gdat(params_equilibrium.shot,'psi_edge','machine',gdat_params.machine); -global_quantities_desc.psi_boundary = 'psi_edge'; -global_quantities.q_95 = gdat(params_equilibrium.shot,'q95','machine',gdat_params.machine); -global_quantities_desc.q_95 = 'q95'; -global_quantities.q_axis = gdat(params_equilibrium.shot,'q0','machine',gdat_params.machine); % will be checked with q_rho? -global_quantities_desc.q_axis = 'q0'; -temp.q_rho = gdat(params_equilibrium.shot,'q_rho','machine',gdat_params.machine); -temp_desc.q_rho = 'q_rho'; +params_eff.data_request = 'li'; +global_quantities.li_3 = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.li_3 = params_eff.data_request; +params_eff.data_request = 'r_axis'; +temp.r_magnetic_axis = gdat(params_equilibrium.shot,params_eff); +temp_desc.r_magnetic_axis = params_eff.data_request; +params_eff.data_request = 'z_axis'; +temp.z_magnetic_axis = gdat(params_equilibrium.shot,params_eff); +temp_desc.z_magnetic_axis = params_eff.data_request; +params_eff.data_request = 'psi_axis'; +temp.psi_axis = gdat(params_equilibrium.shot,params_eff); % needs to add psi_edge sincepsi_axis liuqe assuming 0 dege value +temp_desc.psi_axis = params_eff.data_request; +params_eff.data_request = 'psi_edge'; +global_quantities.psi_boundary = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.psi_boundary = params_eff.data_request; +params_eff.data_request = 'q95'; +global_quantities.q_95 = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.q_95 = params_eff.data_request; +params_eff.data_request = 'q0'; +global_quantities.q_axis = gdat(params_equilibrium.shot,params_eff); % will be checked with q_rho? +global_quantities_desc.q_axis = params_eff.data_request; +params_eff.data_request = 'q_rho'; +temp.q_rho = gdat(params_equilibrium.shot,params_eff); +temp_desc.q_rho = params_eff.data_request; % surface = gdat(params_equilibrium.shot,'surface','machine',gdat_params.machine); % to be added -global_quantities.volume = gdat(params_equilibrium.shot,'volume','machine',gdat_params.machine); -global_quantities_desc.volume = 'volume'; -global_quantities.w_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',gdat_params.machine); -global_quantities_desc.w_mhd = 'w_mhd'; +params_eff.data_request = 'volume'; +global_quantities.volume = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.volume = params_eff.data_request; +params_eff.data_request = 'w_mhd'; +global_quantities.w_mhd = gdat(params_equilibrium.shot,params_eff); +global_quantities_desc.w_mhd = params_eff.data_request; global_quantities_fieldnames = fieldnames(global_quantities); special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments @@ -109,37 +134,49 @@ end % for boundary in addition to lcfs % active_limiter_point = gdat(params_equilibrium.shot,'active_limiter_point','machine',gdat_params.machine); -boundary.elongation = gdat(params_equilibrium.shot,'kappa','machine',gdat_params.machine); -boundary_desc.elongation = 'kappa'; +params_eff.data_request = 'kappa'; +boundary.elongation = gdat(params_equilibrium.shot,params_eff); +boundary_desc.elongation = params_eff.data_request; % elongation_lower = gdat(params_equilibrium.shot,'elongation_lower','machine',gdat_params.machine); % elongation_upper = gdat(params_equilibrium.shot,'elongation_upper','machine',gdat_params.machine); -boundary.minor_radius = gdat(params_equilibrium.shot,'a_minor','machine',gdat_params.machine); -boundary_desc.minor_radius = 'a_minor'; +params_eff.data_request = 'a_minor'; +boundary.minor_radius = gdat(params_equilibrium.shot,params_eff); +boundary_desc.minor_radius = params_eff.data_request; % squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner','machine',gdat_params.machine); % squareness_lower_outer = gdat(params_equilibrium.shot,'squareness_lower_outer','machine',gdat_params.machine); % squareness_upper_inner = gdat(params_equilibrium.shot,'squareness_upper_inner','machine',gdat_params.machine); % squareness_upper_outer = gdat(params_equilibrium.shot,'squareness_upper_outer','machine',gdat_params.machine); % strike_point = gdat(params_equilibrium.shot,'strike_point','machine',gdat_params.machine); -boundary.triangularity = gdat(params_equilibrium.shot,'delta','machine',gdat_params.machine); -boundary_desc.triangularity = 'delta'; -boundary.triangularity_lower = gdat(params_equilibrium.shot,'delta_bottom','machine',gdat_params.machine); -boundary_desc.triangularity_lower = 'delta_bottom'; -boundary.triangularity_upper = gdat(params_equilibrium.shot,'delta_top','machine',gdat_params.machine); -boundary_desc.triangularity_upper = 'delta_top'; -temp.n_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''n_xpts'''',''''liuqe.m'''')','machine',gdat_params.machine); -temp_desc.n_x_point = '''tcv_eq(''''n_xpts'''',''''liuqe.m'''')'''; -temp.r_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''r_xpts'''',''''liuqe.m'''')','machine',gdat_params.machine); -temp_desc.r_x_point = '''tcv_eq(''''r_xpts'''',''''liuqe.m'''')'''; -temp.z_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''z_xpts'''',''''liuqe.m'''')','machine',gdat_params.machine); -temp_desc.z_x_point = '''tcv_eq(''''z_xpts'''',''''liuqe.m'''')'''; -temp.rgeom = gdat(params_equilibrium.shot,'rgeom','machine',gdat_params.machine); -temp_desc.rgeom = 'rgeom'; -temp.zgeom = gdat(params_equilibrium.shot,'zgeom','machine',gdat_params.machine); -temp_desc.zgeom = 'zgeom'; -temp.r_lcfs = gdat(params_equilibrium.shot,'r_contour_edge','machine',gdat_params.machine); -temp_desc.r_lcfs = 'r_contour_edge'; -temp.z_lcfs = gdat(params_equilibrium.shot,'z_contour_edge','machine',gdat_params.machine); -temp_desc.z_lcfs = 'z_contour_edge'; +params_eff.data_request = 'delta'; +boundary.triangularity = gdat(params_equilibrium.shot,params_eff); +boundary_desc.triangularity = params_eff.data_request; +params_eff.data_request = 'delta_bottom'; +boundary.triangularity_lower = gdat(params_equilibrium.shot,params_eff); +boundary_desc.triangularity_lower = params_eff.data_request; +params_eff.data_request = 'delta_top'; +boundary.triangularity_upper = gdat(params_equilibrium.shot,params_eff); +boundary_desc.triangularity_upper = params_eff.data_request; +params_eff.data_request = 'tcv_eq(''''n_xpts'''',''''liuqe.m'''')'; +temp.n_x_point = gdat(params_equilibrium.shot,params_eff); +temp_desc.n_x_point = params_eff.data_request; +params_eff.data_request = 'tcv_eq(''''r_xpts'''',''''liuqe.m'''')'; +temp.r_x_point = gdat(params_equilibrium.shot,params_eff); +temp_desc.r_x_point = params_eff.data_request; +params_eff.data_request = 'tcv_eq(''''z_xpts'''',''''liuqe.m'''')'; +temp.z_x_point = gdat(params_equilibrium.shot,params_eff); +temp_desc.z_x_point = params_eff.data_request; +params_eff.data_request = 'rgeom'; +temp.rgeom = gdat(params_equilibrium.shot,params_eff); +temp_desc.rgeom = params_eff.data_request; +params_eff.data_request = 'zgeom'; +temp.zgeom = gdat(params_equilibrium.shot,params_eff); +temp_desc.zgeom = params_eff.data_request; +params_eff.data_request = 'r_contour_edge'; +temp.r_lcfs = gdat(params_equilibrium.shot,params_eff); +temp_desc.r_lcfs = params_eff.data_request; +params_eff.data_request = 'z_contour_edge'; +temp.z_lcfs = gdat(params_equilibrium.shot,params_eff); +temp_desc.z_lcfs = params_eff.data_request; boundary_fieldnames = fieldnames(boundary); special_fields = {'lcfs', 'geometric_axis', 'x_point'}; % fields needing non-automatic treatments @@ -188,13 +225,19 @@ end % b_min = gdat(params_equilibrium.shot,'b_min','machine',gdat_params.machine); % darea_dpsi = gdat(params_equilibrium.shot,'darea_dpsi','machine',gdat_params.machine); % darea_drho_tor = gdat(params_equilibrium.shot,'darea_drho_tor','machine',gdat_params.machine); -profiles_1d.dpressure_dpsi = gdat(params_equilibrium.shot,'pprime','machine',gdat_params.machine); +params_eff.data_request = 'pprime'; +profiles_1d.dpressure_dpsi = gdat(params_equilibrium.shot,params_eff); +profiles_1d_desc.dpressure_dpsi = params_eff.data_request; % dpsi_drho_tor = gdat(params_equilibrium.shot,'dpsi_drho_tor','machine',gdat_params.machine); % dvolume_dpsi = gdat(params_equilibrium.shot,'dvolume_dpsi','machine',gdat_params.machine); % dvolume_drho_tor = gdat(params_equilibrium.shot,'dvolume_drho_tor','machine',gdat_params.machine); % elongation = gdat(params_equilibrium.shot,'elongation','machine',gdat_params.machine); -profiles_1d.f_df_dpsi = gdat(params_equilibrium.shot,'ttprime','machine',gdat_params.machine); -profiles_1d.f = gdat(params_equilibrium.shot,'rbphi_rho','machine',gdat_params.machine); +params_eff.data_request = 'ttprime'; +profiles_1d.f_df_dpsi = gdat(params_equilibrium.shot,params_eff); +profiles_1d_desc.f_df_dpsi = [params_eff.data_request '* 0.996^2']; +params_eff.data_request = 'rbphi_rho'; +profiles_1d.f = gdat(params_equilibrium.shot,params_eff); +profiles_1d_desc.f = [params_eff.data_request '* 0.996']; profiles_1d.f.data = 0.996 * profiles_1d.f.data; profiles_1d.f_df_dpsi.data = 0.996.^2 .* profiles_1d.f_df_dpsi.data; % geometric_axis = gdat(params_equilibrium.shot,'geometric_axis','machine',gdat_params.machine); @@ -211,14 +254,24 @@ profiles_1d.f_df_dpsi.data = 0.996.^2 .* profiles_1d.f_df_dpsi.data; % j_tor = gdat(params_equilibrium.shot,'j_tor','machine',gdat_params.machine); % magnetic_shear = gdat(params_equilibrium.shot,'magnetic_shear','machine',gdat_params.machine); % mass_density = gdat(params_equilibrium.shot,'mass_density','machine',gdat_params.machine); -profiles_1d.phi = gdat(params_equilibrium.shot,'phi_tor','machine',gdat_params.machine); +params_eff.data_request = 'phi_tor'; +profiles_1d.phi = gdat(params_equilibrium.shot,params_eff); profiles_1d.phi.data = 0.996 * profiles_1d.phi.data; -profiles_1d.pressure = gdat(params_equilibrium.shot,'pressure','machine',gdat_params.machine); +profiles_1d_desc.phi = [params_eff.data_request '* 0.996']; +params_eff.data_request = 'pressure'; +profiles_1d.pressure = gdat(params_equilibrium.shot,params_eff); +profiles_1d_desc.pressure = params_eff.data_request; % psi = gdat(params_equilibrium.shot,'psi_rho','machine',gdat_params.machine); % (could take from .x of any like rhotor and psi_axis, psi_edge from global_quantities) -profiles_1d.q = gdat(params_equilibrium.shot,'q_rho','machine',gdat_params.machine); -profiles_1d.rho_tor = gdat(params_equilibrium.shot,'rhotor','machine',gdat_params.machine); +params_eff.data_request = 'q_rho'; +profiles_1d.q = gdat(params_equilibrium.shot,params_eff); +profiles_1d_desc.q = params_eff.data_request; +params_eff.data_request = 'rhotor'; +profiles_1d.rho_tor = gdat(params_equilibrium.shot,params_eff); +profiles_1d_desc.rho_tor = params_eff.data_request; %rho_tor_norm = gdat(params_equilibrium.shot,'rhotor_norm','machine',gdat_params.machine); % from rho_tor -profiles_1d.rho_volume_norm = gdat(params_equilibrium.shot,'rhovol','machine',gdat_params.machine); +params_eff.data_request = 'rhovol'; +profiles_1d.rho_volume_norm = gdat(params_equilibrium.shot,params_eff); +profiles_1d_desc.rho_volume_norm = params_eff.data_request; % r_inboard = gdat(params_equilibrium.shot,'r_inboard','machine',gdat_params.machine); % r_outboard = gdat(params_equilibrium.shot,'r_outboard','machine',gdat_params.machine); % squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner','machine',gdat_params.machine); @@ -229,7 +282,9 @@ profiles_1d.rho_volume_norm = gdat(params_equilibrium.shot,'rhovol','machine',gd % trapped_fraction = gdat(params_equilibrium.shot,'trapped_fraction','machine',gdat_params.machine); % triangularity_lower = gdat(params_equilibrium.shot,'triangularity_lower','machine',gdat_params.machine); % triangularity_upper = gdat(params_equilibrium.shot,'triangularity_upper','machine',gdat_params.machine); -profiles_1d.volume = gdat(params_equilibrium.shot,'volume_rho','machine',gdat_params.machine,'liuqe',gdat_params.liuqe); +params_eff.data_request = 'volume_rho'; +profiles_1d.volume = gdat(params_equilibrium.shot,params_eff); +profiles_1d_desc.volume = params_eff.data_request; liuqe_opt = gdat_params.liuqe; % default at this stage but could take from gdat params like error bar switch liuqe_opt @@ -240,7 +295,7 @@ switch liuqe_opt case {12,13}, psitbx_str=['LIUQE' num2str(mod(liuqe_opt,10))]; otherwise, error(['Unknown LIUQE version, liuqe = ' num2str(liuqe_opt)]); end -fsd = psitbxtcv2(shot,profiles_1d.volume.t,'FS',psitbx_str); +fsd = psitbxtcv2(shot,profiles_1d.volume.t,'FS',psitbx_str); % will get automatically the correct time interval grho_metric_3D = metric(fsd,-1); % Introduced new anonymous function to compute FS average ... metric_FS = metric(grho_metric_3D.grid,[2,3]); @@ -249,10 +304,12 @@ FS_av = @(x) sum(x.*metric_FS./grho_metric_3D,[2,3])./denom; R=metric(fsd,3); Rm2av=FS_av(1./R.^2); profiles_1d.gm1.data = Rm2av.x; +profiles_1d_desc.gm1 = ['psitbxtcv2 with ' psitbx_str ' then Rm2av=FS_av(1./R.^2)']; %tmp_gm = FS_av(grho_metric_3D.^2./R.^2); % this gives (grad rhopol/R)^2 not gm2 which is grad rhotor^2 %profiles_1d.gm2.data = tmp_gm.x; tmp_gm = FS_av(1./R.^1); profiles_1d.gm9.data = tmp_gm.x; +profiles_1d_desc.gm9 = 'FS_av(1./R.^1)'; tmp_gm = FS_av(grho_metric_3D.^2./R.^2); % grad rhopol^2 to get <grad psi^2> for it=1:numel(ids_equilibrium.time) @@ -322,7 +379,9 @@ profiles_2d.grid_type.description = 'Cylindrical R,Z ala eqdsk'; % j_parallel = gdat(params_equilibrium.shot,'j_parallel','machine',gdat_params.machine); % j_tor = gdat(params_equilibrium.shot,'j_tor','machine',gdat_params.machine); % phi = gdat(params_equilibrium.shot,'phi','machine',gdat_params.machine); -profiles_2d.psi = gdat(params_equilibrium.shot,'psi','machine',gdat_params.machine); % add psi_bound in a second step in special cases +params_eff.data_request = 'psi'; +profiles_2d.psi = gdat(params_equilibrium.shot,params_eff); % add psi_bound in a second step in special cases +profiles_2d_desc.psi = [params_eff.data_request ' adding psi_bound in a 2nd step']; % r = gdat(params_equilibrium.shot,'r','machine',gdat_params.machine); % not to be filled since in grid.dim1 % theta = gdat(params_equilibrium.shot,'theta','machine',gdat_params.machine); % z = gdat(params_equilibrium.shot,'z','machine',gdat_params.machine); % not to be filled since in grid.dim2 diff --git a/matlab/TCV_IMAS/tcv_get_ids_magnetics.m b/matlab/TCV_IMAS/tcv_get_ids_magnetics.m index 2e302cce..e1491454 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_magnetics.m +++ b/matlab/TCV_IMAS/tcv_get_ids_magnetics.m @@ -10,6 +10,8 @@ if exist('gdat_params') 'gdat_params',gdat_params,varargin{:}); else [ids_magnetics, params_magnetics] = tcv_ids_headpart(shot, ids_magnetics_empty,'magnetics','homogeneous_time',0,varargin{:}); + aa=gdat_tcv; + gdat_params = aa.gdat_params; % to get default params end [ids_magnetics.bpol_probe,ids_magnetics_description.bpol_probe]= tcv_ids_bpol_probe(params_magnetics.shot, ids_magnetics.bpol_probe(1),gdat_params); diff --git a/matlab/TCV_IMAS/tcv_get_ids_nbi.m b/matlab/TCV_IMAS/tcv_get_ids_nbi.m index 3b833944..1c2ca51b 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_nbi.m +++ b/matlab/TCV_IMAS/tcv_get_ids_nbi.m @@ -10,7 +10,11 @@ if exist('gdat_params') [ids_nbi, params_nbi] = tcv_ids_headpart(shot,ids_nbi_empty,'nbi','homogeneous_time',0,'gdat_params',gdat_params,varargin{:}); else [ids_nbi, params_nbi] = tcv_ids_headpart(shot,ids_nbi_empty,'nbi','homogeneous_time',0,varargin{:}); + aa=gdat_tcv; + gdat_params = aa.gdat_params; % to get default params end +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) ids_nbi_description=''; % As a general rule, for a new substructure under the main ids, construct a local structure like: @@ -59,26 +63,30 @@ beamlets_group.position(1).phi = 211.9535*pi/180.; beamlets_group.position(1).r = 4.5889; beamlets_group.position(1).z = 0.; beamlets_group.position(2).phi = 295.2416*pi/180.; -beamlets_group.position(2).r = 4.9274; +beamlets_group.position(2).r = 4.9274; beamlets_group.position(2).z = 0.; +params_eff = params_eff_ref; for iunit=1:nb_units ids_nbi.unit{iunit}.identifier = unit_identifier{iunit}; ids_nbi.unit{iunit}.name = unit_name{iunit}; %% power - pow=gdat_tcv(shot,['\results::' results_subname{iunit} ':powr_tcv']); + params_eff.data_request = ['\results::' results_subname{iunit} ':powr_tcv']; + pow=gdat_tcv(shot,params_eff); ids_nbi.unit{iunit}.power_launched.data = pow.data*1e6; ids_nbi.unit{iunit}.power_launched.time = pow.t; - ids_nbi_description.unit{iunit}.power_launched = ['from \results::' results_subname{iunit} ':powr_tcv']; + ids_nbi_description.unit{iunit}.power_launched = params_eff.data_request; %% energy - en=gdat_tcv(shot,['\results::' results_subname{iunit} ':energy']); + params_eff.data_request = ['\results::' results_subname{iunit} ':energy']; + en=gdat_tcv(shot,params_eff); ids_nbi.unit{iunit}.energy.data = en.data*1e3; ids_nbi.unit{iunit}.energy.time = en.t; - ids_nbi_description.unit{iunit}.energy = ['from \results::' results_subname{iunit} ':energy']; + ids_nbi_description.unit{iunit}.energy = params_eff.data_request; %% power & current fractions - p_frac=gdat(shot,['\results::' results_subname{iunit} ':fraction']); + params_eff.data_request = ['\results::' results_subname{iunit} ':fraction']; + p_frac=gdat(shot,params_eff); ids_nbi.unit{iunit}.beam_power_fraction.time = p_frac.t; - ids_nbi_description.unit{iunit}.beam_power_fraction = ['from \results::' results_subname{iunit} ':fraction']; + ids_nbi_description.unit{iunit}.beam_power_fraction = params_eff.data_request; if ~isempty(p_frac.data) && size(p_frac.data,2)>=3 ids_nbi.unit{iunit}.beam_power_fraction.data = p_frac.data(:,1:3)'*0.01; i_frac = p_frac.data(:,1:3).*repmat([1 2 3],size(p_frac.data,1),1); % to be compatible with older matlab version .*[1 2 3] not ok diff --git a/matlab/TCV_IMAS/tcv_get_ids_pf_active.m b/matlab/TCV_IMAS/tcv_get_ids_pf_active.m index e893ee03..a8978b73 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_pf_active.m +++ b/matlab/TCV_IMAS/tcv_get_ids_pf_active.m @@ -8,6 +8,8 @@ if exist('gdat_params') [ids_pf_active, params] = tcv_ids_headpart(shot, ids_pf_active_empty,'pf_active','homogeneous_time',0,'gdat_params',gdat_params,varargin{:}); else [ids_pf_active, params] = tcv_ids_headpart(shot, ids_pf_active_empty,'pf_active','homogeneous_time',0,varargin{:}); + aa=gdat_tcv; + gdat_params = aa.gdat_params; % to get default params end % Get subfield diff --git a/matlab/TCV_IMAS/tcv_get_ids_tf.m b/matlab/TCV_IMAS/tcv_get_ids_tf.m index d658a82a..4ce27ae1 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_tf.m +++ b/matlab/TCV_IMAS/tcv_get_ids_tf.m @@ -13,11 +13,14 @@ ids_tf = ids_tf_empty; if exist('gdat_params') [ids_tf, params_tf] = tcv_ids_headpart(shot, ids_tf_empty,'tf','homogeneous_time',0,'gdat_params',gdat_params,varargin{:}); + params_eff = gdat_params; + params_eff.data_request = 'b0'; + tmp = gdat_tcv(shot,params_eff); else [ids_tf, params_tf] = tcv_ids_headpart(shot, ids_tf_empty,'tf','homogeneous_time',0,varargin{:}); + tmp = gdat_tcv(shot, 'b0'); end -tmp = gdat_tcv(shot, 'b0'); if ~ischar(tmp.data) ids_tf.r0 = tmp.r0; ids_tf_description.r0 = ' b0.r0 from gdat_tcv(shot, ''b0'')'; diff --git a/matlab/TCV_IMAS/tcv_get_ids_wall.m b/matlab/TCV_IMAS/tcv_get_ids_wall.m index 9d90d36e..19404fea 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_wall.m +++ b/matlab/TCV_IMAS/tcv_get_ids_wall.m @@ -9,9 +9,11 @@ if exist('gdat_params') [ids_wall, params] = tcv_ids_headpart(shot, ids_wall_empty,'wall','homogeneous_time',0,'gdat_params',gdat_params,varargin{:}); else [ids_wall, params] = tcv_ids_headpart(shot, ids_wall_empty,'wall','homogeneous_time',0,varargin{:}); + aa=gdat_tcv; + gdat_params = aa.gdat_params; % to get default params end -% Get data of outline +% Get data of outline, only static, no need for gdat_params and time_out parameter for example [ids_wall.description_2d,ids_wall_description.description_2d]= tcv_ids_wall_description_2d(params.shot, ids_wall.description_2d(1)); % make arrays not filled in empty: not the case for magnetics diff --git a/matlab/TCV_IMAS/tcv_ids_bpol_probe.m b/matlab/TCV_IMAS/tcv_ids_bpol_probe.m index 5027b8e8..d5764521 100644 --- a/matlab/TCV_IMAS/tcv_ids_bpol_probe.m +++ b/matlab/TCV_IMAS/tcv_ids_bpol_probe.m @@ -14,13 +14,20 @@ error_bar = 'delta'; if exist('gdat_params') && isfield(gdat_params,'error_bar') && ~isempty(gdat_params.error_bar) error_bar = gdat_params.error_bar; end +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) % Get data -tmp = gdat_tcv( shot, '\MAGNETICS::BPOL_003'); +params_eff = params_eff_ref; +params_eff.data_request = '\MAGNETICS::BPOL_003'; +tmp = gdat_tcv( shot,params_eff); +tmp_desc = params_eff.data_request; names = tmp.dim{2}; time = tmp.dim{1}; data = tmp.data; -ang_m = gdat_tcv( shot, 'static("ANG_M")'); +params_eff.data_request = 'static("ANG_M")'; +ang_m = gdat_tcv(shot,params_eff); +ang_m_desc = params_eff.data_request; % from mapping file from Holger: probes_name={ '001' '002' '003' '004' '005' '006' '007' '008' '009' '010' '011' '012' '013' '014' '015' '016' '017' '018' '019' ... @@ -52,7 +59,8 @@ for ii=1:Nprobes sub_ids_struct_out{ii}.field.time = time; sub_ids_struct_out{ii}.toroidal_angle = 0.; % to see if should match sector 3 (bpol003) - aa = gdat_tcv( shot, ['static("ANG_M")[$1]'',''' sub_ids_struct_out{ii}.name '']); + params_eff.data_request = ['static("ANG_M")[$1]'',''' sub_ids_struct_out{ii}.name '']; + aa = gdat_tcv(shot,params_eff); sub_ids_struct_out{ii}.poloidal_angle = -aa.data; ij=strmatch(names{ii},probes_name); sub_ids_struct_out{ii}.area = probes_area(ij); diff --git a/matlab/TCV_IMAS/tcv_ids_circuit.m b/matlab/TCV_IMAS/tcv_ids_circuit.m index 3b83cd4c..6009c46d 100644 --- a/matlab/TCV_IMAS/tcv_ids_circuit.m +++ b/matlab/TCV_IMAS/tcv_ids_circuit.m @@ -14,12 +14,15 @@ error_bar = 'delta'; if exist('gdat_params') && isfield(gdat_params,'error_bar') && ~isempty(gdat_params.error_bar) error_bar = gdat_params.error_bar; end +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) %% Get power supply/coils names for each circuit. [tcv_circuit_info] = tcv_ids_pf_active_definition(); % Preallocate memory and get data ids_struct_out(1:tcv_circuit_info.ntotcircuits) = ids_structures; +params_eff = params_eff_ref; for ii=1:tcv_circuit_info.ntotcircuits if shot == -1 % model shot % replace by dummy @@ -28,7 +31,8 @@ for ii=1:tcv_circuit_info.ntotcircuits warning('no time data loaded for shot %d',shot); ids_struct_out_description{ii}.current = 'not loaded'; else - tmpdata = gdat_tcv(shot,['' tcv_circuit_info.mds_paths{ii} '']); % Get current + params_eff.data_request = ['' tcv_circuit_info.mds_paths{ii} '']; + tmpdata = gdat_tcv(shot,params_eff);; ids_struct_out_description{ii}.current = ['from ' tmpdata.data_fullpath]; end ids_struct_out{ii}.current.data = tmpdata.data; diff --git a/matlab/TCV_IMAS/tcv_ids_coil.m b/matlab/TCV_IMAS/tcv_ids_coil.m index 1ef3f60f..13eb86f0 100644 --- a/matlab/TCV_IMAS/tcv_ids_coil.m +++ b/matlab/TCV_IMAS/tcv_ids_coil.m @@ -16,6 +16,8 @@ error_bar = 'delta'; if exist('gdat_params') && isfield(gdat_params,'error_bar') && ~isempty(gdat_params.error_bar) error_bar = gdat_params.error_bar; end +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) % Coils that can be characterized by R, Z and a % crosssectional area are described as distinct coils with a single element, this also corresponds to coils with distinct TCV names). @@ -27,7 +29,6 @@ end % dealt with in the circuit description and in the machine mapping. % In practice the dissipated energy in the pf ' % coils is not a relevant limit. -mdsopen(shot) %% Get power supply/coils names for each circuit. [pf_def] = tcv_ids_pf_active_definition(); @@ -38,6 +39,7 @@ ncircuits2ids = numel(coil_names2ids); %% Get geometrical data from the static tree % Information from the static tree - Z-zeroed with respect to coils +mdsopen(shot) r_c = mdsvalue('_r_c=static(''r_c'')'); % R position z_c = mdsvalue('_z_c=static(''z_c'')'); % Z position w_c = mdsvalue('static(''w_c'')'); % width rectangular description @@ -65,6 +67,7 @@ h_c(iG)=sqrt(xsect_c(iG)); %% Put data to ids structure ids_struct_out(1:sizepfc) = ids_structures; ind_coil_ids = 0; +params_eff = params_eff_ref; for ii=1:ncircuits2ids ncoil2ids = numel(coil_names2ids{ii}); % number of coils for a given circuit @@ -72,22 +75,26 @@ for ii=1:ncircuits2ids ind_coil_ids = ind_coil_ids +1; ids_struct_out{ind_coil_ids}.name = coil_names2ids{ii}{jj}; ids_struct_out_description{ind_coil_ids}.name = ['through aa=tcv_ids_pf_active_definition from aa.coil_names ']; - + % time-varying data if shot == -1 tmpdata.dim{1} = []; tmpdata.data = []; warning('no time data loaded for shot %d',shot); else - tmpdata = tdi(mds_paths2ids{ii}); + params_eff.data_request = mds_paths2ids{ii}; + tmpdata = gdat_tcv(shot,params_eff); end % $$$ ids_struct_out{ind_coil_ids}.current.data = tmpdata.data; % $$$ ids_struct_out_description{ind_coil_ids}.current.data = ['from ' mds_paths2ids{ii}]; ids_struct_out{ind_coil_ids}.current.data = pf_def.coil_current_signs{ii}(jj) .* tmpdata.data; ids_struct_out_description{ind_coil_ids}.current.data = ['from ' num2str(pf_def.coil_current_signs{ii}(jj)) ' * ' mds_paths2ids{ii}]; +try ids_struct_out{ind_coil_ids}.current.time = tmpdata.dim{1}; +catch + keyboard + end - % Find index on static tree tmpind = find(strcmp(ids_struct_out{ind_coil_ids}.name, namepfc)); ids_struct_out{ind_coil_ids}.element{1}.geometry.rectangle.r = r_c(tmpind); @@ -102,7 +109,7 @@ for ii=1:ncircuits2ids ids_struct_out_description{ind_coil_ids}.element{1}.turns_with_sign = ['from static(''nt_c'')']; ids_struct_out{ind_coil_ids}.element{1}.geometry.geometry_type = 2; % 1 outline, 2 rectangle, 4 arc of circle ids_struct_out_description{ind_coil_ids}.element{1}.geometry.geometry_type = 'rectangle'; - + ids_struct_out{ind_coil_ids}.resistance = res_c(tmpind); ids_struct_out_description{ind_coil_ids}.resistance = ['from static(''res_c'')']; end diff --git a/matlab/TCV_IMAS/tcv_ids_flux_loop.m b/matlab/TCV_IMAS/tcv_ids_flux_loop.m index f5fe9e9d..5ba7c87c 100644 --- a/matlab/TCV_IMAS/tcv_ids_flux_loop.m +++ b/matlab/TCV_IMAS/tcv_ids_flux_loop.m @@ -14,9 +14,14 @@ error_bar = 'delta'; if exist('gdat_params') && isfield(gdat_params,'error_bar') && ~isempty(gdat_params.error_bar) error_bar = gdat_params.error_bar; end +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) % Get data -tmp = gdat_tcv(shot, 'tcv_idealloop("FLUX")'); +params_eff = params_eff_ref; +params_eff.data_request = 'tcv_idealloop("FLUX")'; +tmp = gdat_tcv(shot,params_eff); +tmp_desc = params_eff.data_request; names = tmp.dim{2}; time = tmp.dim{1}; data = tmp.data; diff --git a/matlab/TCV_IMAS/tcv_ids_headpart.m b/matlab/TCV_IMAS/tcv_ids_headpart.m index 911f2baf..d1d45e6e 100644 --- a/matlab/TCV_IMAS/tcv_ids_headpart.m +++ b/matlab/TCV_IMAS/tcv_ids_headpart.m @@ -92,13 +92,11 @@ if ~isempty(params_ids_generic.gdat_params) ids_generic.ids_properties.comment(end+1:end+length(comment_add)+2) = ['; ' comment_add]; end else - warning('cocos_in and cocos_out expected to be fields of gdat_params') + disp('cocos_in and cocos_out expected to be fields of gdat_params') end if isfield(params_ids_generic.gdat_params,'error_bar') provider_add = ['; error_bar option: ' params_ids_generic.gdat_params.error_bar]; ids_generic.ids_properties.provider(end+1:end+length(provider_add)) = provider_add; - else - warning('cocos_in and cocos_out expected to be fields of gdat_params') end end ids_generic.ids_properties.homogeneous_time = params_ids_generic.homogeneous_time; diff --git a/matlab/TCV_IMAS/tcv_ids_ip.m b/matlab/TCV_IMAS/tcv_ids_ip.m index e843cb0a..fe1e786e 100644 --- a/matlab/TCV_IMAS/tcv_ids_ip.m +++ b/matlab/TCV_IMAS/tcv_ids_ip.m @@ -14,12 +14,19 @@ error_bar = 'delta'; if exist('gdat_params') && isfield(gdat_params,'error_bar') && ~isempty(gdat_params.error_bar) error_bar = gdat_params.error_bar; end +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) % Get data -tmp = gdat_tcv( shot, 'ip_trapeze'); +params_eff = params_eff_ref; +params_eff.data_request = 'ip_trapeze'; +tmp = gdat_tcv(shot,params_eff); +tmp_desc = params_eff.data_request; time = tmp.dim{1}; data = tmp.data; -tmpdml = gdat_tcv( shot, '\results::dmlcor'); +params_eff.data_request = '\results::dmlcor'; +tmpdml = gdat_tcv(shot,params_eff); +tmpdml_desc = params_eff.data_request; % Preallocate dimension ids_struct_out = ids_structures; diff --git a/matlab/TCV_IMAS/tcv_ids_supply.m b/matlab/TCV_IMAS/tcv_ids_supply.m index 0776a12a..6fd1b29e 100644 --- a/matlab/TCV_IMAS/tcv_ids_supply.m +++ b/matlab/TCV_IMAS/tcv_ids_supply.m @@ -14,12 +14,15 @@ error_bar = 'delta'; if exist('gdat_params') && isfield(gdat_params,'error_bar') && ~isempty(gdat_params.error_bar) error_bar = gdat_params.error_bar; end +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) %% Get power supply/coils names for each circuit. [tcv_circuit_info] = tcv_ids_pf_active_definition(); % Preallocate memory and get data ids_struct_out(1:tcv_circuit_info.ntotpowersupplies) = ids_structures; +params_eff = params_eff_ref; for ii=1:tcv_circuit_info.ntotpowersupplies if shot ==-1 warning('TCV:IDS:NoTimeData','no time data loaded for shot %d',shot); @@ -29,7 +32,8 @@ for ii=1:tcv_circuit_info.ntotpowersupplies ids_struct_out_description{ii}.current = sprintf('Not loaded for shot %d',shot); else - tmpdata = gdat_tcv(shot,['' tcv_circuit_info.mds_paths{ii} '']); % Get current + params_eff.data_request = ['' tcv_circuit_info.mds_paths{ii} '']; + tmpdata = gdat_tcv(shot,params_eff);; ids_struct_out_description{ii}.current = ['abs value from ' tmpdata.data_fullpath]; end diff --git a/matlab/gdat2time_out.m b/matlab/gdat2time_out.m new file mode 100644 index 00000000..16e4f442 --- /dev/null +++ b/matlab/gdat2time_out.m @@ -0,0 +1,288 @@ +function [gdat_data_out] = gdat2time_out(gdat_data_in,option,others_as_data,varargin); +% +% function [gdat_data_out] = gdat2time_out(gdat_data_in,time_out,option,gdat_timedim,varargin); +% +% Aim: Reduce data output to desired time interval or compute values at desired time_instances +% +% Inputs: +% gdat_data_in: gdat_data structure +% option: 1: standard .data, .t, .dim{gdat_timedim} to be changed +% 21, 22: deal with grids_1d as well: 21 with rhopol 1D, 22 with rhopol 2D +% +% others_as_data: {'field1','field2',...} other fields gdat_data_in.(others_as_data{i}) to be treated as .data +% +% Inferred/necessary inputs: +% gdat_data_in.gdat_params.time_out: time interval [t1,t2] or list of times to get output +% gdat_data_in.mapping_for.(gdat_data_in.gdat_params.machine).gdat_timedim: dim of time +% +% Outputs: +% gdat_data_out: with reduced data +% + +gdat_data_out = gdat_data_in; + +if isempty(gdat_data_in) || isempty(gdat_data_in.t) || ~isnumeric(gdat_data_in.data) || ~isnumeric(gdat_data_in.t) + return +end + +if ~exist('option') || isempty(option) + option = 1; +end + +iothers_as_data = false; +others_as_data_eff = ''; +if exist('others_as_data') && ~isempty(others_as_data) && iscell(others_as_data) + others_as_data_eff = others_as_data; + iothers_as_data = true; +end + +try + dim_time = gdat_data_in.mapping_for.(gdat_data_in.gdat_params.machine).gdat_timedim; +catch ME + try + getReport(ME) + error(['expect gdat_timedim in mapping_for.' gdat_data_in.gdat_params.machine ' from data_request: ' gdat_data_in.gdat_request]); + catch ME2 + getReport(ME2) + error(['expect gdat_timedim in mapping_for.machine fro data_request']); + end +end + +try + time_out = sort(gdat_data_in.gdat_params.time_out); +catch ME + getReport(ME) + error(['expect gdat_data_in.gdat_params.time_out']) +end + +gdat_data_out = gdat_data_in; + +if any(option == [1, 21, 22]) + if numel(time_out) == 2 + % time interval provided, take all existing time points within this interval + ij = [gdat_data_out.t>=time_out(1) & gdat_data_out.t<=time_out(2)]; + gdat_data_out.t = gdat_data_out.t(ij); + if isfield(gdat_data_out,'dim'); gdat_data_out.dim{dim_time} = gdat_data_out.t; end + nbdims = numel(size(gdat_data_out.data)); + nbdims_eff = nbdims; + if nbdims==2 && any(size(gdat_data_out.data)==1) + nbdims_eff = 1; + end + switch nbdims_eff + case 1 + gdat_data_out.data = gdat_data_out.data(ij); + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(ij);']); + end + end + case 2 + if dim_time ==1 + gdat_data_out.data = gdat_data_out.data(ij,:); + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']); + if ~iscell(atmp) + eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(ij,:);']); + else + eval(['gdat_data_out.' others_as_data_eff{j} '{1} = gdat_data_out.' others_as_data_eff{j} '{1}(ij);']); + end + end + end + else + gdat_data_out.data = gdat_data_out.data(:,ij); + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']); + if ~iscell(atmp) + eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(:,ij);']); + else + eval(['gdat_data_out.' others_as_data_eff{j} '{2} = gdat_data_out.' others_as_data_eff{j} '{2}(ij);']); + end + end + end + end + case 3 + if dim_time == 1 + gdat_data_out.data = gdat_data_out.data(ij,:,:); + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']); + if ~iscell(atmp) + eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(ij,:,:);']); + else + eval(['gdat_data_out.' others_as_data_eff{j} '{1} = gdat_data_out.' others_as_data_eff{j} '{1}(ij);']); + end + end + end + elseif dim_time == 2 + gdat_data_out.data = gdat_data_out.data(:,ij,:); + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']); + if ~iscell(atmp) + eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(:,ij,:);']); + else + eval(['gdat_data_out.' others_as_data_eff{j} '{2} = gdat_data_out.' others_as_data_eff{j} '{2}(ij);']); + end + end + end + else + gdat_data_out.data = gdat_data_out.data(:,:,ij); + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']); + if ~iscell(atmp) + eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_out.' others_as_data_eff{j} '(:,:,ij);']); + else + eval(['gdat_data_out.' others_as_data_eff{j} '{3} = gdat_data_out.' others_as_data_eff{j} '{3}(ij);']); + end + end + end + end + otherwise + warning(['nbdims_eff = ' num2str(nbdims_eff) ' not yet defined in gdat2time_out, ask O. Sauter if needed']); + end + else + % Assume 1 or several (>2) time points are requested. Interpolate to get data at these time instances specifically + % default time_out_interp='linear', will implement {'interpos','tension','option',..} later when needed + % time interval provided, take all existing time points within this interval + ij = iround_os(gdat_data_out.t,time_out); + gdat_data_out.t = time_out; + if isfield(gdat_data_out,'dim'); gdat_data_out.dim{dim_time} = gdat_data_out.t; end + nbdims = numel(size(gdat_data_out.data)); + nbdims_eff = nbdims; + if nbdims==2 && any(size(gdat_data_out.data)==1) + nbdims_eff = 1; + end + switch nbdims_eff + case 1 + gdat_data_out.data = interp1(gdat_data_in.t,gdat_data_in.data,gdat_data_out.t); + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['gdat_data_out.' others_as_data_eff{j} ' = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} ',gdat_data_out.t);']); + end + end + case 2 + len_t = numel(ij); + if dim_time ==1 + gdat_data_out.data = NaN*ones(len_t,size(gdat_data_out.data,2)); + for ix=1:size(gdat_data_out.data,2) + gdat_data_out.data([1:len_t],ix) = interp1(gdat_data_in.t,gdat_data_in.data(:,ix),gdat_data_out.t); + end + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(len_t,size(gdat_data_out.' others_as_data_eff{j} ',2));']); + end + for j=1:length(others_as_data_eff) + ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']); + for ix=1:ix_len + eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(:,ix),gdat_data_out.t);']); + end + end + end + else + gdat_data_out.data = NaN*ones(size(gdat_data_out.data,1),len_t); + for ix=1:size(gdat_data_out.data,1) + gdat_data_out.data(ix,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.data(ix,:),gdat_data_out.t); + end + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(size(gdat_data_out.' others_as_data_eff{j} ',1),len_t);']); + end + for j=1:length(others_as_data_eff) + ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',1);']); + for ix=1:ix_len + eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(ix,:),gdat_data_out.t);']); + end + end + end + end + case 3 + len_t = numel(ij); + if dim_time == 1 + gdat_data_out.data = NaN*ones(len_t,size(gdat_data_out.data,2),size(gdat_data_out.data,3)); + for ix=1:size(gdat_data_out.data,2) + for jy=1:size(gdat_data_out.data,3) + gdat_data_out.data([1:len_t],ix,jy) = interp1(gdat_data_in.t,gdat_data_in.data(:,ix,jy),gdat_data_out.t); + end + end + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(len_t,size(gdat_data_out.' others_as_data_eff{j} ... + ',2),size(gdat_data_out.' others_as_data_eff{j} ',3));']); + end + for j=1:length(others_as_data_eff) + ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']); + jy_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',3);']); + for ix=1:ix_len + for jy=1:jy_len + eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix,jy) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(:,ix,jy),gdat_data_out.t);']); + end + end + end + end + elseif dim_time == 2 + gdat_data_out.data = NaN*ones(size(gdat_data_out.data,1),len_t,size(gdat_data_out.data,3)); + for ix=1:size(gdat_data_out.data,1) + for jy=1:size(gdat_data_out.data,3) + gdat_data_out.data(ix,[1:len_t],jy) = interp1(gdat_data_in.t,gdat_data_in.data(ix,:,jy),gdat_data_out.t); + end + end + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(size(gdat_data_out.' others_as_data_eff{j} ... + ',1),len_t,size(gdat_data_out.' others_as_data_eff{j} ',3));']); + end + for j=1:length(others_as_data_eff) + ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',1);']); + jy_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',3);']); + for ix=1:ix_len + for jy=1:jy_len + eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t],jy) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(ix,:,jy),gdat_data_out.t);']); + end + end + end + end + else + gdat_data_out.data = NaN*ones(size(gdat_data_out.data,1),size(gdat_data_out.data,2),len_t); + for ix=1:size(gdat_data_out.data,1) + for jy=1:size(gdat_data_out.data,2) + gdat_data_out.data(ix,jy,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.data(ix,jy,:),gdat_data_out.t); + end + end + if iothers_as_data + for j=1:length(others_as_data_eff) + eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(size(gdat_data_out.' others_as_data_eff{j} ... + ',1),size(gdat_data_out.' others_as_data_eff{j} ',2),len_t);']); + end + for j=1:length(others_as_data_eff) + ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',1);']); + jy_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']); + for ix=1:ix_len + for jy=1:jy_len + eval(['gdat_data_out.' others_as_data_eff{j} '(ix,jy,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(ix,jy,:),gdat_data_out.t);']); + end + end + end + end + end + otherwise + warning(['nbdims_eff = ' num2str(nbdims_eff) ' not yet defined in gdat2time_out, ask O. Sauter if needed']); + end + end +else + error(['option = ' num2str(option) ' not yet implemented']) +end + +if any(option==[21,22]) + % add extra fields then correct + ab_tmp_rho = gdat_data_out.grids_1d; + ab_tmp_rho.data = gdat_data_out.grids_1d.rhotornorm; ab_tmp_rho.t = gdat_data_out.t; + ab_tmp_rho.gdat_params = gdat_data_out.gdat_params; ab_tmp_rho.mapping_for = gdat_data_out.mapping_for; + extra_fields = {'psi','rhotornorm','rhovolnorm','rhopolnorm_equil','rhotornorm_equil', ... + 'rhovolnorm_equil','rhotor_edge','volume_edge','b0'}; + if option==22; extra_fields{end+1} = 'rhopolnorm'; end + [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,extra_fields); + gdat_data_out.grids_1d = rmfield(ab_tmp_rho,{'data','t','gdat_params','mapping_for'}); +end -- GitLab