From 5fe3f50b86407e26e0b54debbb60747b9153b69e Mon Sep 17 00:00:00 2001 From: Olivier Sauter <olivier.sauter@epfl.ch> Date: Tue, 2 May 2017 06:05:51 +0000 Subject: [PATCH] start adding D3D, halpha and {'spectroscopy','\fs04'} works git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@7291 d63d8f72-b253-0410-a779-e742ad2e26cf --- crpptbx/D3D/d3d_help_parameters.m | 71 + crpptbx/D3D/d3d_requests_mapping.m | 387 ++++ crpptbx/D3D/gdat_d3d.m | 1971 +++++++++++++++++ crpptbx/D3D/get_signal_d3d.m | 28 + crpptbx/gdat.m | 2 + crpptbx/gdat_plot.m | 2 +- crpptbx/gdatpaths.m | 1 + .../get_data_request_names_from_gdat_xxx.m | 2 +- 8 files changed, 2462 insertions(+), 2 deletions(-) create mode 100644 crpptbx/D3D/d3d_help_parameters.m create mode 100644 crpptbx/D3D/d3d_requests_mapping.m create mode 100644 crpptbx/D3D/gdat_d3d.m create mode 100644 crpptbx/D3D/get_signal_d3d.m diff --git a/crpptbx/D3D/d3d_help_parameters.m b/crpptbx/D3D/d3d_help_parameters.m new file mode 100644 index 00000000..ec8d86a7 --- /dev/null +++ b/crpptbx/D3D/d3d_help_parameters.m @@ -0,0 +1,71 @@ +function help_struct = d3d_help_parameters(parameter_list) +% +% retrieve from present table the relevant help lines for the parameter_list{:} +% should come from sqlite database at some point... +% +% return the whole help structure if parameter_list empty or not provided +% +% do: +% help_struct = d3d_help_parameters(fieldnames(gdat_data.gdat_params)); +% +% to get relevant help description +% + +% Defaults +help_struct_all = struct(... + 'data_request', ['automatically filled in by gdat, name of request used in gdat call.' char(10) ... + 'contains current list of keywords if gdat called with no arguments: aa=gdat;' char(10) ... + 'Note shot value should not be in params so params can be used to load same data from another shot'] ... + ,'machine', 'machine name like ''TCV'', ''D3D'', case insensitive' ... + ,'doplot', '0 (default), if 1 calls gdat_plot for a new figure, -1 plot over current figure with hold all, see gdat_plot for details' ... + ,'liuqe','liuqe version 1 (default), 2, 3 for LIUQE1, 2, 3 resp. or -1 for model values' ... + ,'nverbose','1 (default) displays warnings, 0: only errors, >=3: displays all extra information' ... + ); + +% D3D 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' ... +% $$$ ' 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, nete: 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.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']; +% $$$ help_struct_all.edge = '0 (default), 1 to get edge Thomson values'; +help_struct_all.fit = '0, no fit profiles, 1 (default) if fit profiles desired as well, relevant for _rho profiles. See also fit_type'; +help_struct_all.fit_type = 'type of fits ''std'' (default) uses diagnostic error bars, ''pedestal'', uses manual error bars with smaller values outside 0.8'; +help_struct_all.fit_nb_rho_points = 'nb of points for the radial mesh over which the fits are evaluated for the fitted profiles, it uses an equidistant mesh at this stage'; +help_struct_all.source = 'sxr: ''G'' (default, with ssx), camera name ''J'', ''G'', ...[[F-M], case insensitive;; cxrs: ''CEZ'' (default), ''CMZ'''; +help_struct_all.camera = ['[] (default, all), [i1 i2 ...] chord nbs ([1 3 5] if only chords 1, 3 and 5 are desired), ''central'' uses J_049']; +help_struct_all.freq = '''slow'', default (means ssx, 500kHz), lower sampling; ''fast'' full samples 2MHz; integer value nn for downsampling every nn''th points'; +%help_struct_all. = ''; + +%help_struct_all. = ''; + + + +if ~exist('parameter_list') || isempty(parameter_list) + help_struct = help_struct_all; + return +end + +if iscell(parameter_list) + parameter_list_eff = parameter_list; +elseif ischar(parameter_list) + % assume only one parameter provided as char string + parameter_list_eff{1} = parameter_list; +else + disp(['unknown type of parameter_list in d3d_help_parameters.m']) + parameter_list + return +end + +fieldnames_all = fieldnames(help_struct_all); +for i=1:length(parameter_list_eff) + if any(strcmp(fieldnames_all,parameter_list_eff{i})) + help_struct.(parameter_list_eff{i}) = help_struct_all.(parameter_list_eff{i}); + end +end diff --git a/crpptbx/D3D/d3d_requests_mapping.m b/crpptbx/D3D/d3d_requests_mapping.m new file mode 100644 index 00000000..71e14c86 --- /dev/null +++ b/crpptbx/D3D/d3d_requests_mapping.m @@ -0,0 +1,387 @@ +function mapping = d3d_requests_mapping(data_request) +% +% Information pre-defined for gdat_d3d, you can add cases here to match official cases in gdat_d3d, allowing backward compatibility +% + +% Defaults +mapping = struct(... + 'label', '', ... + 'method', '', ... + 'expression','', ... + 'timedim', -1, ... % dim which is the time is the database, to copy in .t, the other dims are in .x (-1 means last dimension) + 'gdat_timedim',[], ... % if need to reshape data and dim orders to have timedim as gdat_timedim (shifting time to gdat_timedim) + 'min', -inf, ... + 'max', inf); +% Note that gdat_timedim is set to timedim at end of this function if empty +% gdat_timedim should have effective index of the time dimension in gdat + +if ~exist('data_request') || isempty(data_request) + return +end + +% default label: data_request keyword itself +mapping.label = data_request; + +% for D3D, following choices are set so far: +% method = 'signal' then expression contains the shotfile, diagnostic and if needed the experiment +% expression is a cell array +% method = 'expression', then expression is the expression to return gdat_tmp... +% method = 'switchcase', then there will be a specific case within gdat_d3d (usual case when not directly a signal) +% +% label is used for plotting +if iscell(data_request) % || (~ischar(data_request) && length(data_request)>1) + mapping.label = data_request; + mapping.method = 'signal'; % assume a full tracename is given, so just try with tdi (could check there are some ":", etc...) + mapping.expression = data_request; + mapping.gdat_timedim = mapping.timedim; + return +end + +switch lower(data_request) + case 'a_minor' + mapping.timedim = 1; + mapping.label = 'a\_minor'; + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''r_inboard'';' ... + 'gdat_tmp=gdat_d3d(shot,params_eff);gdat_tmp.r_inboard=gdat_tmp.data;' ... + 'params_eff.data_request=''r_outboard'';' ... + 'gdat_tmp2=gdat_d3d(shot,params_eff);gdat_tmp.r_outboard=gdat_tmp2.data;' ... + 'gdat_tmp.data = 0.5.*(gdat_tmp2.data-gdat_tmp.data);gdat_tmp.label=''' mapping.label ''';' ... + 'gdat_tmp.gdat_request=''' data_request ''';']; + case 'b0' + mapping.timedim = 1; + mapping.label = 'B_0'; + mapping.method = 'signal'; + mapping.expression = [{'FPC'},{'BTF'}]; + case 'beta' + mapping.timedim = 1; + mapping.label = '\beta'; + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''betan'';' ... + 'gdat_tmp=gdat_d3d(shot,params_eff);' ... + 'params_eff.data_request=''ip'';gdat_tmp2=gdat_d3d(shot,params_eff);' ... + 'params_eff.data_request=''b0'';gdat_tmp3=gdat_d3d(shot,params_eff);' ... + 'params_eff.data_request=''a_minor'';gdat_tmp4=gdat_d3d(shot,params_eff);' ... + 'tmp_data_ip=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ... + 'tmp_data_b0=interp1(gdat_tmp3.t,gdat_tmp3.data,gdat_tmp.t,[],NaN);' ... + 'tmp_data_a=interp1(gdat_tmp4.t,gdat_tmp4.data,gdat_tmp.t,[],NaN);' ... + 'gdat_tmp.data = 0.01.*abs(gdat_tmp.data.*tmp_data_ip./1e6./tmp_data_a./tmp_data_b0);']; + case 'betan' + mapping.timedim = 1; + mapping.label = '\beta_N'; + mapping.method = 'signal'; + mapping.expression = [{'TOT'},{'beta_N'}]; + % in many cases, in particular just after an experiment, betaN is not present in TOT, thus compute it from 2/3Wmhd/V /(B0^2/2mu0) + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=[{''TOT''},{''beta_2N''}];' ... + 'gdat_tmp=gdat_d3d(shot,params_eff); if isempty(gdat_tmp.data);' ... + 'params_eff.data_request=''ip'';gdat_ip=gdat_d3d(shot,params_eff);' ... + 'params_eff.data_request=''b0'';gdat_b0=gdat_d3d(shot,params_eff);' ... + 'params_eff.data_request=''a_minor'';gdat_aminor=gdat_d3d(shot,params_eff);' ... + 'params_eff.data_request=''wmhd'';gdat_tmp=gdat_d3d(shot,params_eff);' ... + 'params_eff.data_request=''volume'';gdat_vol=gdat_d3d(shot,params_eff);' ... + 'tmp_data_ip=interp1(gdat_ip.t,gdat_ip.data,gdat_tmp.t,[],NaN);' ... + 'tmp_data_b0=interp1(gdat_b0.t,gdat_b0.data,gdat_tmp.t,[],NaN);' ... + 'tmp_data_a=interp1(gdat_aminor.t,gdat_aminor.data,gdat_tmp.t,[],NaN);' ... + 'tmp_data_vol=interp1(gdat_vol.t,gdat_vol.data,gdat_tmp.t,[],NaN);' ... + 'gdat_tmp.data = 100.*abs(2./3.*gdat_tmp.data./tmp_data_vol.*8e-7.*pi./tmp_data_b0.^2./tmp_data_ip.*1e6.*tmp_data_a.*tmp_data_b0);end;']; + case 'betap' + mapping.timedim = 1; + mapping.label = '\beta_p'; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'betpol'}]; + case {'cxrs', 'cxrs_rho'} + mapping.timedim = 2; + mapping.label = 'cxrs'; + mapping.method = 'switchcase'; + mapping.expression = ''; + case 'delta' + mapping.timedim = 1; + mapping.label = 'delta'; + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''delta_bottom''; ' ... + 'gdat_tmp=gdat_d3d(shot,params_eff);params_eff.data_request=''delta_top'';' ... + 'gdat_tmp2=gdat_d3d(shot,params_eff);gdat_tmp.data = 0.5.*(gdat_tmp.data+gdat_tmp2.data);']; + case 'delta_top' + mapping.label = 'delta\_top'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'delRoben'}]; + case 'delta_bottom' + mapping.label = 'delta\_bottom'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'delRuntn'}]; + case {'ece', 'eced', 'ece_rho', 'eced_rho'} + mapping.timedim = 2; + mapping.method = 'switchcase'; + mapping.expression = ''; + case 'eqdsk' + mapping.timedim = 2; + mapping.method = 'switchcase'; % could use function make_eqdsk directly? + mapping.expression = ''; + case 'equil' + mapping.gdat_timedim = 2; + mapping.method = 'switchcase'; % could use function make_eqdsk directly? + mapping.expression = ''; + case 'halpha' + mapping.timedim = 1; + mapping.label = 'Halpha'; + mapping.method = 'signal'; + mapping.expression = [{'SPECTROSCOPY'},{'\fs04'}]; + case 'h_scalings' + mapping.label = 'H_{scal}'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'TTH'},{'H/L-facs'}]; + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=[{''TTH''},{''Wmhd''}];' ... + 'gdat_tmp=gdat_d3d(shot,params_eff);ij=find(~isnan(gdat_tmp.data)); gdat_tmp.data_raw=gdat_tmp.data;' ... + 'tmp_data=interpos(gdat_tmp.t(ij),gdat_tmp.data(ij),gdat_tmp.t,-1e5);' ... + 'gdat_tmp.data = max(tmp_data,0.);']; + case 'ioh' + mapping.timedim = 1; + mapping.label = 'I ohmic transformer'; + mapping.method = 'signal'; + mapping.expression = [{'MBI'},{'IOH'}]; + case 'ip' + mapping.timedim = 1; + mapping.label = 'Plasma current'; + mapping.method = 'signal'; + mapping.expression = [{'MAG'},{'Ipa'}]; + case 'kappa' + mapping.timedim = 1; + mapping.label = '\kappa'; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'k'}]; + case 'kappa_top' + mapping.timedim = 1; + mapping.label = '\kappa^{top}'; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'koben'}]; + case 'kappa_bottom' + mapping.timedim = 1; + mapping.label = '\kappa_{bottom}'; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'kuntn'}]; + case 'li' + mapping.timedim = 1; + mapping.label = 'l_i'; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'li'}]; + case 'mhd' + mapping.timedim = 1; + mapping.label = 'Odd and Even n'; + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request={''MOD'',''OddN''}; ' ... + 'gdat_tmp=gdat_d3d(shot,params_eff);gdat_tmp.data=reshape(gdat_tmp.data,length(gdat_tmp.data),1 );' ... + 'gdat_tmp.dim{1}=gdat_tmp.t;gdat_tmp.dim{2}=[1 2];gdat_tmp.x=gdat_tmp.dim{2};' ... + 'gdat_tmp.n_odd.data = gdat_tmp.data;gdat_tmp.n_odd.data_request=params_eff.data_request;' ... + 'params_eff.data_request={''MOD'',''EvenN''};' ... + 'gdat_tmp2=gdat_d3d(shot,params_eff);gdat_tmp.data(:,2)=reshape(gdat_tmp2.data,length(gdat_tmp2.data),1);' ... + 'gdat_tmp.n_even.data = gdat_tmp2.data;gdat_tmp.n_even.data_request=params_eff.data_request;gdat_tmp.label={''n\_odd'',''n\_even''};' ... + 'params_eff.data_request={''MOD'',''OddNAmp''};' ... + 'gdat_tmp2=gdat_d3d(shot,params_eff);gdat_tmp.n_odd.amp=reshape(gdat_tmp2.data,length(gdat_tmp2.data),1);' ... + 'gdat_tmp.n_odd.amp_t=gdat_tmp2.t;' ... + 'params_eff.data_request={''MOD'',''EvenNAmp''};' ... + 'gdat_tmp2=gdat_d3d(shot,params_eff);gdat_tmp.n_even.amp=reshape(gdat_tmp2.data,length(gdat_tmp2.data),1);' ... + 'gdat_tmp.n_even.amp_t=gdat_tmp2.t;' ... + 'gdat_tmp.full_path=''MOD/Odd in data and .n_odd; .n_even'';' ... + 'gdat_tmp.gdat_request=''mhd'';gdat_tmp.gdat_params.data_request=gdat_tmp.gdat_request;']; + case 'ne' + mapping.timedim = 2; + mapping.method = 'switchcase'; + case 'neint' + mapping.timedim = 1; + mapping.label = 'line integrated el. density'; + mapping.method = 'signal'; + mapping.expression = [{'DCN'},{'H-1'}]; + case 'nel' + mapping.timedim = 1; + mapping.label = 'line-averaged el. density'; + mapping.expression = [{'FPG'},{'lenH-1'}]; + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''neint'';' ... + 'gdat_tmp=gdat_d3d(shot,params_eff);params_eff.data_request=[{''FPG''},{''lenH-1''}];' ... + 'gdat_tmp2=gdat_d3d(shot,params_eff);ij=find(gdat_tmp2.data==0);gdat_tmp2.data(ij)=NaN;' ... + 'tmp_data=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ... + 'gdat_tmp.data = gdat_tmp.data./(tmp_data+1e-5);']; + case 'ne_rho' + mapping.timedim = 2; + mapping.label = 'ne'; + mapping.method = 'switchcase'; + case 'nete_rho' + mapping.timedim = 2; + mapping.label = 'ne and Te'; + mapping.method = 'switchcase'; + case 'ni' + mapping.method = 'switchcase'; % especially since might have option fit, etc + case 'pgyro' + mapping.timedim = 1; + mapping.label = 'EC gyros'; + mapping.method = 'switchcase'; + case 'powers' + mapping.timedim = 1; + mapping.label = 'various powers'; + mapping.method = 'switchcase'; + case 'psi_axis' + mapping.timedim = 1; + mapping.method = 'switchcase'; % there is psi_axis-psi_edge in FPG but otherwise complicated to get from equil, thus needs swticth case + mapping.label ='psi_\axis' ; + case 'psi_edge' + mapping.timedim = 1; + mapping.method = 'switchcase'; % is set to zero, so not in tree nodes + mapping.label = 'psi\_edge'; + case 'q0' + mapping.timedim = 1; + mapping.label = 'q_0'; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'q0'}]; + case 'q95' + mapping.timedim = 1; + mapping.label = 'q_{95}'; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'q95'}]; + case 'q_edge' + mapping.timedim = 1; + mapping.label = 'q_{edge}}'; + mapping.method = 'expression'; + mapping.method = 'switchcase'; + mapping.expression = [{'FPG'},{'q95'}]; + mapping.expression = []; + case 'q_rho' + mapping.timedim = 2; + mapping.gdat_timedim = 2; + mapping.label = 'q'; + mapping.method = 'switchcase'; + case 'rgeom' + mapping.label = 'Rgeom'; + mapping.timedim = 1; + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''r_inboard'';' ... + 'gdat_tmp=gdat_d3d(shot,params_eff);gdat_tmp.r_inboard=gdat_tmp.data;' ... + 'params_eff.data_request=''r_outboard'';' ... + 'gdat_tmp2=gdat_d3d(shot,params_eff);gdat_tmp.r_outboard=gdat_tmp2.data;' ... + 'gdat_tmp.data = 0.5.*(gdat_tmp2.data+gdat_tmp.data);gdat_tmp.label=''' mapping.label ''';' ... + 'gdat_tmp.gdat_request=''' data_request ''';']; + case 'r_inboard' + mapping.label = 'R\_inboard'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'Rin'}]; + case 'r_outboard' + mapping.label = 'R\_outboard'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'Raus'},{'D3DD'}]; + case 'rhotor' + mapping.timedim = 2; + mapping.method = 'switchcase'; + mapping.label = 'rhotor'; + case 'rhotor_edge' + mapping.timedim = 1; + mapping.method = 'switchcase'; + mapping.label = 'rhotor\_edge'; + case 'rhotor_norm' + mapping.timedim = 1; + mapping.method = 'switchcase'; + mapping.label = 'rhotor\_norm'; + case 'rhovol' + mapping.timedim = 2; + mapping.label = 'rhovol\_norm'; + mapping.method = 'switchcase'; + case 'rmag' + mapping.label = 'R\_magaxis'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'Rmag'},{'D3DD'}]; + case 'sxr' + mapping.timedim = 1; + mapping.gdat_timedim = 2; + mapping.method = 'switchcase'; + case 'tau_tot' + mapping.label = '\tau_{tot}'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'TOT'},{'tau_tot'},{'D3DD'}]; + case 'te' + mapping.timedim = 2; + mapping.label = 'Te'; + mapping.method = 'switchcase'; + case 'te_rho' + mapping.timedim = 2; + mapping.label = 'Te'; + mapping.method = 'switchcase'; + case 'ti' + mapping.label = 'Ti'; + mapping.method = 'switchcase'; + case 'vloop' + mapping.label = 'Vloop'; + mapping.timedim = 1; + % mapping.method = 'signal'; + % mapping.expression = [{'MAG'},{'ULid12'},{'D3DD'}]; + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=[{''MAG''},{''ULid12''},{''D3DD''}];' ... + 'gdat_tmp=gdat_d3d(shot,params_eff);ij=find(~isnan(gdat_tmp.data));' ... + 'tmp_data=interpos(gdat_tmp.t,gdat_tmp.data,-3e4);' ... + 'gdat_tmp.data_smooth = tmp_data;gdat_tmp.gdat_request=''vloop'';gdat_tmp.gdat_params.data_request=gdat_tmp.gdat_request;']; + case 'volume' + mapping.label = 'Volume'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'Vol'},{'D3DD'}]; + case 'volume_rho' + mapping.timedim = 2; + mapping.label = 'volume\_norm'; + mapping.method = 'switchcase'; + case 'wmhd' + mapping.label = 'Wmhd'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'Wmhd'},{'D3DD'}]; + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=[{''FPG''},{''Wmhd''},{''D3DD''}];' ... + 'gdat_tmp=gdat_d3d(shot,params_eff);ij=find(~isnan(gdat_tmp.data)); gdat_tmp.data_raw=gdat_tmp.data;' ... + 'tmp_data=interpos(gdat_tmp.t(ij),gdat_tmp.data(ij),gdat_tmp.t,-1e5);' ... + 'gdat_tmp.data = max(tmp_data,0.);']; + case 'zeff' + mapping.label = 'zeff from cxrs'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'ZES'},{'Zeff'},{'D3DD'}]; + case 'zgeom' + mapping.label = 'Zgeom'; + mapping.timedim = 1; + mapping.method = 'expression'; + mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=[{''FPG''},{''Zoben''},{''D3DD''}];' ... + 'gdat_tmp=gdat_d3d(shot,params_eff);gdat_tmp.z_top=gdat_tmp.data;' ... + 'params_eff.data_request=[{''FPG''},{''Zunt''},{''D3DD''}];' ... + 'gdat_tmp2=gdat_d3d(shot,params_eff);gdat_tmp.z_bottom=gdat_tmp2.data;' ... + 'gdat_tmp.data = 0.5.*(gdat_tmp2.data+gdat_tmp.data);gdat_tmp.label=''' mapping.label ''';' ... + 'gdat_tmp.gdat_request=''' data_request ''';']; + + case 'zmag' + mapping.label = 'Z\_magaxis'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'FPG'},{'Zmag'},{'D3DD'}]; + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % extra D3D cases (not necessarily in official data_request name list) + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + case 'transp' + mapping.label = 'transp output'; + mapping.method = 'switchcase'; + + + otherwise + mapping.label = data_request; + mapping.method = 'signal'; % assume a full tracename is given, so just try with tdi (could check there are some ":", etc...) + mapping.expression = data_request; + +end + +if isempty(mapping.gdat_timedim) + mapping.gdat_timedim = mapping.timedim; +end diff --git a/crpptbx/D3D/gdat_d3d.m b/crpptbx/D3D/gdat_d3d.m new file mode 100644 index 00000000..1b640907 --- /dev/null +++ b/crpptbx/D3D/gdat_d3d.m @@ -0,0 +1,1971 @@ +function [gdat_data,gdat_params,error_status,varargout] = gdat_d3d(shot,data_request,varargin) +% +% function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin) +% +% Aim: get data from a given machine using full path or keywords. +% data_request are and should be case independent (transformed in lower case in the function and outputs) +% +% If no inputs are provided, return the list of available pre-defined data_request in gdat_data and default parameters gdat_params +% +% Inputs: +% +% no inputs: return default parameters in a structure form in gdat_params +% shot: shot number +% data_request: keyword (like 'ip') or trace name or structure containing all parameters but shot number +% varargin{i},varargin{i+1},i=1:nargin-2: additional optional parameters given in pairs: param_name, param_value +% The optional parameters list might depend on the data_request +% examples of optional parameters: +% 'plot',1 (plot is set by default to 0) +% 'machine','D3D' (the default machine is the local machine) +% +% +% Outputs: +% +% gdat_data: structure containing the data, the time trace if known and other useful information +% gdat_data.t : time trace +% gdat_data.data: requested data values +% gdat_data.dim : values of the various coordinates related to the dimensions of .data(:,:,...) +% note that one of the dim is the time, replicated in .t for clarity +% gdat_data.dimunits : units of the various dimensions, 'dimensionless' if dimensionless +% gdat_data.error_bar : if provided +% gdat_data.gdat_call : list of parameters provided in the gdat call (so can be reproduced) +% gdat_data.shot: shot number +% gdat_data.machine: machine providing the data +% gdat_data.gdat_request: keyword for gdat if relevant +% gdat_data.data_fullpath: full path to the data node if known and relevant, or relevant expression called if relevant +% gdat_data.gdat_params: copy gdat_params for completeness (gdat_params contains a .help structure detailing each parameter) +% gdat_data.xxx: any other relevant information +% +% gdat_params contains the options relevant for the called data_request. It also contains a help structure for each option +% eg.: param1 in gdat_params.param1 and help in gdat_params.help.param1 +% +% Examples: +% (should add working examples for various machines (provides working shot numbers for each machine...)) +% +% [a1,a2]=gdat; +% a2.data_request = 'Ip'; +% a3=gdat(32827,a2); % gives input parameters as a structure, allows to call the same for many shots +% a4=gdat('opt1',123,'opt2',[1 2 3],'shot',48832,'data_request','Ip','opt3','aAdB'); % all in pairs +% a5=gdat(32827,'ip'); % standard call +% a6=gdat(32827,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (note all lowercase in output) + +% +% Comments for local developer: +% This gdat is just a "header routine" calling the gdat for the specific machine gdat_`machine`.m which can be called +% directly, thus which should be able to treat the same type of input arguments +% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Prepare some variables, etc + +varargout{1}=cell(1,1); +error_status=1; + +% construct main default parameters structure +gdat_params.data_request = ''; +default_machine = 'd3d'; + +gdat_params.machine=default_machine; +gdat_params.doplot = 0; +gdat_params.exp_name = 'D3DD'; +gdat_params.nverbose = 1; + +% construct list of keywords from global set of keywords and specific D3D set +% get data_request names from centralized table +%%% data_request_names = get_data_request_names; % do not use xlsx anymore but scan respective machine_requests_mapping.m files +data_request_names = get_data_request_names_from_gdat_xxx(default_machine); +% add D3D specific to all: +if ~isempty(data_request_names.d3d) + d3d_names = fieldnames(data_request_names.d3d); + for i=1:length(d3d_names) + data_request_names.all.(d3d_names{i}) = data_request_names.d3d.(d3d_names{i}); + end +end +data_request_names_all = sort(fieldnames(data_request_names.all)); + +% construct default output structure +gdat_data.data = []; +gdat_data.units = []; +gdat_data.dim = []; +gdat_data.dimunits = []; +gdat_data.t = []; +gdat_data.x = []; +gdat_data.shot = []; +gdat_data.gdat_request = []; +gdat_data.gdat_params = gdat_params; +gdat_data.data_fullpath = []; + + +% Treat inputs: +ivarargin_first_char = 3; +data_request_eff = ''; +if nargin>=2 && ischar(data_request); data_request = lower(data_request); end + +gdat_data.gdat_request = data_request_names_all; % so if return early gives list of possible request names +gdat_data.gdat_params.help = d3d_help_parameters(fieldnames(gdat_data.gdat_params)); +gdat_params = gdat_data.gdat_params; + +% no inputs +if nargin==0 + % return defaults and list of keywords + return +end + +do_mdsopen_mdsclose = 1; +% treat 1st arg +if nargin>=1 + if isempty(shot) + % means mdsopen(shot) already performed + shot=-999; % means not defined + do_mdsopen_mdsclose = 0; + elseif isnumeric(shot) + gdat_data.shot = shot; + elseif ischar(shot) + ivarargin_first_char = 1; + else + if gdat_params.nverbose>=1; warning('type of 1st argument unexpected, should be numeric or char'); end + error_status=2; + return + end + if nargin==1 + % Only shot number given. If there is a default data_request set it and continue, otherwise return + return + end +end +% 2nd input argument if not part of pairs +if nargin>=2 && ivarargin_first_char~=1 + if isempty(data_request) + return + end + % 2nd arg can be a structure with all options except shot_number, or a string for the pathname or keyword, or the start of pairs string/value for the parameters + if isstruct(data_request) + if ~isfield(data_request,'data_request') + if gdat_params.nverbose>=1; warning('expects field data_request in input parameters structure'); end + error_status=3; + return + end + %data_request.data_request = lower(data_request.data_request); + data_request_eff = data_request.data_request; + gdat_params = data_request; + else + % since data_request is char check from nb of args if it is data_request or a start of pairs + if mod(nargin-1,2)==0 + ivarargin_first_char = 2; + else + ivarargin_first_char = 3; + data_request_eff = data_request; + end + end +end +if ~isstruct(data_request) + gdat_params.data_request = data_request_eff; +end + +% if start pairs from shot or data_request, shift varargin +if ivarargin_first_char==1 + varargin_eff{1} = shot; + varargin_eff{2} = data_request; + varargin_eff(3:nargin) = varargin(:); +elseif ivarargin_first_char==2 + varargin_eff{1} = data_request; + varargin_eff(2:nargin-1) = varargin(:); +else + varargin_eff(1:nargin-2) = varargin(:); +end + +% extract parameters from pairs of varargin: +if (nargin>=ivarargin_first_char) + if mod(nargin-ivarargin_first_char+1,2)==0 + for i=1:2:nargin-ivarargin_first_char+1 + if ischar(varargin_eff{i}) + % enforce lower case for any character driven input + if ischar(varargin_eff{i+1}) + gdat_params.(lower(varargin_eff{i})) = lower(varargin_eff{i+1}); + else + gdat_params.(lower(varargin_eff{i})) = varargin_eff{i+1}; + end + else + if gdat_params.nverbose>=1; warning(['input argument nb: ' num2str(i) ' is incorrect, expects a character string']); end + error_status=401; + return + end + end + else + if gdat_params.nverbose>=1; warning('number of input arguments incorrect, cannot make pairs of parameters'); end + error_status=402; + return + end +end +data_request_eff = gdat_params.data_request; % in case was defined in pairs +% default equil: +if ~isfield(gdat_params,'equil'); gdat_params.equil = 'EQI'; end +if isfield(gdat_params,'source') && (any(strcmp(lower(gdat_params.source),'ide')) || any(strcmp(lower(gdat_params.source),'idg')) ... + || any(strcmp(lower(gdat_params.source),'ida'))); gdat_params.equil = 'IDE'; end + +% if it is a request_keyword can obtain description: +if ischar(data_request_eff) || length(data_request_eff)==1 + ij=strmatch(data_request_eff,data_request_names_all,'exact'); +else + ij=[]; +end + +if ~isempty(ij); + gdat_data.gdat_request = data_request_names_all{ij}; + if isfield(data_request_names.all.(data_request_names_all{ij}),'description') && ~isempty(data_request_names.all.(data_request_names_all{ij}).description) + % copy description of keyword + gdat_data.request_description = data_request_names.all.(data_request_names_all{ij}).description; + end +else + if ~isempty(data_request_eff); gdat_data.gdat_request = data_request_eff; end +end + +% special treatment if shot and data_request given within pairs +if isfield(gdat_params,'shot') + shot = gdat_params.shot; % should use only gdat_params.shot but change shot to make sure + gdat_data.shot = gdat_params.shot; + gdat_params=rmfield(gdat_params,'shot'); +end + +if ~isfield(gdat_params,'data_request') || isempty(gdat_params.data_request) + % warning('input for ''data_request'' missing from input arguments') % might be correct, asking for list of requests + error_status=5; + return +end +gdat_data.gdat_params = gdat_params; + +% re-assign main variables to make sure use the one in gdat_data structure +shot = gdat_data.shot; +data_request_eff = gdat_data.gdat_params.data_request; +error_status = 6; % at least reached this level + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Specifications on how to get the data provided in d3d_requests_mapping +mapping_for_d3d = d3d_requests_mapping(data_request_eff); +gdat_data.label = mapping_for_d3d.label; + +% fill again at end to have full case, but here to have present status in case of early return +gdat_data.gdat_params.help = d3d_help_parameters(fieldnames(gdat_data.gdat_params)); +gdat_data.mapping_for.d3d = mapping_for_d3d; +gdat_params = gdat_data.gdat_params; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 1st treat the simplest method: "signal" +if strcmp(mapping_for_d3d.method,'signal') + % allow changing main source with simple signals, 'source' parameter relates to mapping_for_d3d.expression{1} + if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || ~ischar(gdat_data.gdat_params.source) + gdat_data.gdat_params.source = mapping_for_d3d.expression{1}; + else + end + [shoteff,stat] = mdsopen(['atlas.gat.com::' mapping_for_d3d.expression{1}],shot); + if shoteff~=shot + disp('cannot open shot'); + return + end + [gdat_data,error_status]=get_signal_d3d(mapping_for_d3d.expression{2},gdat_data); + if mapping_for_d3d.timedim==-1 + mapping_for_d3d.timedim = 1; + mapping_for_d3d.gdat_timedim = 1; + end + if error_status~=0 + if gdat_params.nverbose>=3; disp(['error after get_signal_d3d in signal with data_request_eff= ' data_request_eff]); end + return + end + + gdat_data.data_fullpath=mapping_for_d3d.expression; + % end of method "signal" + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +elseif strcmp(mapping_for_d3d.method,'expression') + % 2nd: method="expression" + % assume expression contains an expression to evaluate and which creates a local structure into variable gdat_tmp + % we copy the structure, to make sure default nodes are defined and to avoid if return is an closed object like tdi + % eval_expr = [mapping_for_d3d.expression ';']; + gdata_data_orig = gdat_data; + fields_to_copy = {'gdat_params', 'gdat_request', 'label', 'data_fullpath', 'request_description', 'mapping_for'}; + eval([mapping_for_d3d.expression ';']); + for i=1:length(fields_to_copy) + if isfield(gdata_data_orig,fields_to_copy{i}) + gdat_tmp.(fields_to_copy{i}) = gdata_data_orig.(fields_to_copy{i}); + end + end + if isempty(gdat_tmp) || (~isstruct(gdat_tmp) & ~isobject(gdat_tmp)) + if (gdat_params.nverbose>=1); warning(['expression does not create a gdat_tmp structure: ' mapping_for_d3d.expression]); end + error_status=801; + return + end + tmp_fieldnames = fieldnames(gdat_tmp); + if sum(strcmp(tmp_fieldnames,'data'))==0 % note: cannot do isfield since gdat_tmp might be an object + if (gdat_params.nverbose>=1); warning(['expression does not return a child name ''data'' for ' data_request_eff]); end + end + for i=1:length(tmp_fieldnames) + gdat_data.(tmp_fieldnames{i}) = gdat_tmp.(tmp_fieldnames{i}); + end + % add .t and .x in case only dim is provided + % do not allow shifting of timedim since should be treated in the relevant expression + ijdim=find(strcmp(tmp_fieldnames,'dim')==1); + if ~isempty(ijdim) + nbdims = length(gdat_data.dim); + if mapping_for_d3d.timedim==-1; + mapping_for_d3d.timedim = nbdims; + if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_d3d.timedim = nbdims-1; end + end + dim_nontim = setdiff([1:nbdims],mapping_for_d3d.timedim); + ijt=find(strcmp(tmp_fieldnames,'t')==1); + if isempty(ijt) + gdat_data.t = gdat_data.dim{mapping_for_d3d.timedim}; + end + ijx=find(strcmp(tmp_fieldnames,'x')==1); + if isempty(ijx) + if ~isempty(dim_nontim) + % since most cases have at most 2d, copy as array if data is 2D and as cell if 3D or more + if length(dim_nontim)==1 + gdat_data.x = gdat_data.dim{dim_nontim(1)}; + else + gdat_data.x = gdat_data.dim(dim_nontim); + end + end + end + gdat_data.data_fullpath=mapping_for_d3d.expression; + end + % end of method "expression" + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +elseif strcmp(mapping_for_d3d.method,'switchcase') + switch data_request_eff % not lower(...) since data_request_eff should be lower case already at this stage + case {'cxrs', 'cxrs_rho'} + % if 'fit' option is added: 'fit',1, then the fitted profiles are returned which means "_rho" is effectively called + if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit) || ~isnumeric(gdat_data.gdat_params.fit) + gdat_data.gdat_params.fit = 0; + elseif gdat_data.gdat_params.fit==1 + data_request_eff = 'cxrs_rho'; + end + exp_name_eff = gdat_data.gdat_params.exp_name; + if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || strcmp(upper(gdat_data.gdat_params.source),'CEZ') + gdat_data.gdat_params.source = 'CEZ'; + r_node = 'R_time'; + z_node = 'z_time'; + elseif strcmp(upper(gdat_data.gdat_params.source),'CMZ') + gdat_data.gdat_params.source = 'CMZ'; + r_node = 'R'; + z_node = 'z'; + else + warning(['source = ' gdat_data.gdat_params.source ' not expected with data_request= ' data_request_eff]); + return + end + diag_name = gdat_data.gdat_params.source; + % R, Z positions of measurements + try + [r_time]=sf2ab(diag_name,shot,r_node,'-exp',exp_name_eff); + catch ME_R_time + % assume no shotfile + getReport(ME_R_time,'basic'); + return + end + gdat_data.r = r_time.value{1}; + inotok=find(gdat_data.r<=0); + gdat_data.r(inotok) = NaN; + [z_time]=sf2ab(diag_name,shot,z_node,'-exp',exp_name_eff); + gdat_data.z = z_time.value{1}; + inotok=find(gdat_data.z<=0); + gdat_data.z(inotok) = NaN; + [time]=sf2tb(diag_name,shot,'time','-exp',exp_name_eff); + gdat_data.t = time.value; + gdat_data.dim{1} = {gdat_data.r , gdat_data.z}; + gdat_data.dimunits{1} = 'R, Z [m]'; + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits{2} = 't [s]'; + gdat_data.x = gdat_data.dim{1}; + % vrot + [a,error_status]=rdaD3D_eff(shot,diag_name,'vrot',exp_name_eff); + if isempty(a.data) || isempty(a.t) || error_status>0 + if gdat_params.nverbose>=3; + a + disp(['with data_request= ' data_request_eff]) + end + end + gdat_data.vrot.data = a.data; + if any(strcmp(fieldnames(a),'units')); gdat_data.vrot.units=a.units; end + if any(strcmp(fieldnames(a),'name')); gdat_data.vrot.name=a.name; end + gdat_data.vrot.label = 'vrot_tor'; + [aerr,e]=rdaD3D_eff(shot,diag_name,'err_vrot',exp_name_eff); + gdat_data.vrot.error_bar = aerr.data; + % Ti + % [a,e]=rdaD3D_eff(shot,diag_name,'Ti',exp_name_eff); + % [aerr,e]=rdaD3D_eff(shot,diag_name,'err_Ti',exp_name_eff); + [a,e]=rdaD3D_eff(shot,diag_name,'Ti_c',exp_name_eff); + [aerr,e]=rdaD3D_eff(shot,diag_name,'err_Ti_c',exp_name_eff); + gdat_data.ti.data = a.data; + gdat_data.data = a.data; + gdat_data.label = [gdat_data.label '/Ti']; + if any(strcmp(fieldnames(a),'units')); gdat_data.ti.units=a.units; end + if any(strcmp(fieldnames(a),'units')); gdat_data.units=a.units; end + if any(strcmp(fieldnames(a),'name')); gdat_data.ti.name=a.name; end + gdat_data.ti.label = 'Ti_c'; + gdat_data.ti.error_bar = aerr.data; + gdat_data.error_bar = aerr.data; + gdat_data.data_fullpath=[exp_name_eff '/' diag_name '/' a.name ';vrot, Ti in data, {r;z} in .x']; + % + if strcmp(data_request_eff,'cxrs_rho') + params_equil = gdat_data.gdat_params; + params_equil.data_request = 'equil'; + [equil,params_equil,error_status] = gdat_d3d(shot,params_equil); + if error_status>0 + if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end + return + end + gdat_data.gdat_params.equil = params_equil.equil; + gdat_data.equil = equil; + inb_chord_cxrs=size(gdat_data.data,1); + inb_time_cxrs=size(gdat_data.data,2); + psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs); + rhopolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs); + rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs); + rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs); + % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)] + time_equil=[min(gdat_data.t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.t(end)+0.1)]; + iok=find(~isnan(gdat_data.r(:,1))); + for itequil=1:length(time_equil)-1 + rr=equil.Rmesh(:,itequil); + zz=equil.Zmesh(:,itequil); + psirz_in = equil.psi2D(:,:,itequil); + it_cxrs_inequil = find(gdat_data.t>time_equil(itequil) & gdat_data.t<=time_equil(itequil+1)); + if ~isempty(it_cxrs_inequil) + if strcmp(upper(gdat_data.gdat_params.source),'CEZ') + rout=gdat_data.r(iok,it_cxrs_inequil); + zout=gdat_data.z(iok,it_cxrs_inequil); + else + rout=gdat_data.r(iok); + zout=gdat_data.z(iok); + end + psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout); + if strcmp(upper(gdat_data.gdat_params.source),'CEZ') + psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil)); + else + psi_out(iok,it_cxrs_inequil) = repmat(psi_at_routzout,[1,length(it_cxrs_inequil)]); + end + rhopolnorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))); + for it_cx=1:length(it_cxrs_inequil) + rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]); + rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]); + end + end + end + gdat_data.psi = psi_out; + gdat_data.rhopolnorm = rhopolnorm_out; + gdat_data.rhotornorm = rhotornorm_out; + gdat_data.rhovolnorm = rhovolnorm_out; + % defaults for fits, so user always gets std structure + gdat_data.fit.rhotornorm = []; % same for both ti and vrot + gdat_data.fit.rhopolnorm = []; + gdat_data.fit.t = []; + gdat_data.fit.ti.data = []; + gdat_data.fit.ti.drhotornorm = []; + gdat_data.fit.vrot.data = []; + gdat_data.fit.vrot.drhotornorm = []; + gdat_data.fit.raw.rhotornorm = []; + gdat_data.fit.raw.ti.data = []; + gdat_data.fit.raw.vrot.data = []; + fit_tension_default = -1; + if isfield(gdat_data.gdat_params,'fit_tension') + fit_tension = gdat_data.gdat_params.fit_tension; + else + fit_tension = fit_tension_default; + end + if ~isstruct(fit_tension) + fit_tension_eff.ti = fit_tension; + fit_tension_eff.vrot = fit_tension; + fit_tension = fit_tension_eff; + else + if ~isfield(fit_tension,'ti'); fit_tension.ti = fit_tension_default; end + if ~isfield(fit_tension,'vrot '); fit_tension.vrot = fit_tension_default; end + end + gdat_data.gdat_params.fit_tension = fit_tension; + if isfield(gdat_data.gdat_params,'fit_nb_rho_points') + fit_nb_rho_points = gdat_data.gdat_params.fit_nb_rho_points; + else + fit_nb_rho_points = 201; + end + gdat_data.gdat_params.fit_nb_rho_points = fit_nb_rho_points; + % + if gdat_data.gdat_params.fit==1 + % add fits + gdat_data.fit.raw.rhotornorm = NaN*ones(size(gdat_data.ti.data)); + gdat_data.fit.raw.ti.data = NaN*ones(size(gdat_data.ti.data)); + gdat_data.fit.raw.vrot.data = NaN*ones(size(gdat_data.vrot.data)); + rhotornormfit = linspace(0,1,fit_nb_rho_points)'; + gdat_data.fit.rhotornorm = rhotornormfit; + gdat_data.fit.t = gdat_data.t; + for it=1:length(gdat_data.t) + % make rhotor->rhopol transformation for each time since equilibrium might have changed + irhook=find(gdat_data.rhotornorm(:,it)>0 & gdat_data.rhotornorm(:,it)<1); % no need for ~isnan + [rhoeff isort]=sort(gdat_data.rhotornorm(irhook,it)); + gdat_data.fit.rhopolnorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.rhopolnorm(irhook(isort),it); 1],rhotornormfit,-0.1,[2 2],[0 1]); + idata = find(gdat_data.ti.data(:,it)>0 & gdat_data.rhotornorm(:,it)<1.01); + if length(idata)>0 + gdat_data.fit.ti.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it); + gdat_data.fit.ti.raw.data(idata,it) = gdat_data.ti.data(idata,it); + gdat_data.fit.ti.raw.error_bar(idata,it) = gdat_data.ti.error_bar(idata,it); + gdat_data.fit.vrot.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it); + gdat_data.fit.vrot.raw.data(idata,it) = gdat_data.vrot.data(idata,it); + gdat_data.fit.vrot.raw.error_bar(idata,it) = gdat_data.vrot.error_bar(idata,it); + [rhoeff,irhoeff] = sort(gdat_data.rhotornorm(idata,it)); + rhoeff = [0; rhoeff]; + % they are some strange low error_bars, so remove these by max(Ti/error_bar)<=10; and changing it to large error bar + tieff = gdat_data.ti.data(idata(irhoeff),it); + ti_err_eff = gdat_data.ti.error_bar(idata(irhoeff),it); + ij=find(tieff./ti_err_eff>10.); + if ~isempty(ij); ti_err_eff(ij) = tieff(ij)./0.1; end + vroteff = gdat_data.vrot.data(idata(irhoeff),it); + vrot_err_eff = gdat_data.vrot.error_bar(idata(irhoeff),it); + ij=find(vroteff./vrot_err_eff>10.); + if ~isempty(ij); vrot_err_eff(ij) = vroteff(ij)./0.1; end + % + tieff = [tieff(1); tieff]; + ti_err_eff = [1e4; ti_err_eff]; + vroteff = [vroteff(1); vroteff]; + vrot_err_eff = [1e5; vrot_err_eff]; + [gdat_data.fit.ti.data(:,it), gdat_data.fit.ti.drhotornorm(:,it)] = interpos(rhoeff,tieff,rhotornormfit,fit_tension.ti,[1 0],[0 0],ti_err_eff); + [gdat_data.fit.vrot.data(:,it), gdat_data.fit.vrot.drhotornorm(:,it)] = interpos(rhoeff,vroteff,rhotornormfit,fit_tension.vrot,[1 0],[0 0],vrot_err_eff); + end + end + end + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'ece', 'eced', 'ece_rho', 'eced_rho'} + nth_points = 13; + if isfield(gdat_data.gdat_params,'nth_points') && ~isempty(gdat_data.gdat_params.nth_points) + nth_points = gdat_data.gdat_params.nth_points; + else + gdat_data.gdat_params.nth_points = nth_points; + end + channels = -1; + if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels) + channels = gdat_data.gdat_params.channels; + end + if nth_points>=10 + match_rz_to_time = 1; + else + match_rz_to_time = 0; + end + if isfield(gdat_data.gdat_params,'match_rz_to_time') && ~isempty(gdat_data.gdat_params.match_rz_to_time) + match_rz_to_time = gdat_data.gdat_params.match_rz_to_time; + else + gdat_data.gdat_params.match_rz_to_time = match_rz_to_time; + end + time_interval = [-Inf +Inf]; + if isfield(gdat_data.gdat_params,'time_interval') && ~isempty(gdat_data.gdat_params.time_interval) + time_interval = gdat_data.gdat_params.time_interval; + else + gdat_data.gdat_params.time_interval = time_interval; + end + % + if isfield(gdat_data.gdat_params,'diag_name') && ~isempty(gdat_data.gdat_params.diag_name) + diag_name = gdat_data.gdat_params.diag_name; + if strcmp(upper(diag_name),'RMD'); gdat_data.gdat_params.exp_name = 'ECED'; end + else + diag_name = 'CEC'; + gdat_data.gdat_params.diag_name = diag_name; + end + exp_name_eff = gdat_data.gdat_params.exp_name; + if strcmp(data_request_eff,'eced') + exp_name_eff = 'ECED'; + gdat_data.gdat_params.exp_name = exp_name_eff; + end + [a,e]=rdaD3D_eff(shot,diag_name,'Trad-A',exp_name_eff,time_interval); + % [a,e]=rdaD3D_eff(shot,diag_name,'Trad-A',exp_name_eff); + inb_chord = size(a.data,1); + if channels(1)<=0 + channels = [1:inb_chord]; + end + gdat_data.dim{1} = channels; + gdat_data.gdat_params.channels = channels; + if nth_points>1 + gdat_data.data = a.data(channels,1:nth_points:end); + gdat_data.dim{2} = a.t(1:nth_points:end); + else + gdat_data.data = a.data(channels,:); + gdat_data.dim{2} = a.t; + end + gdat_data.x = gdat_data.dim{1}; + gdat_data.t = gdat_data.dim{2}; + gdat_data.dimunits=[{'channels'} ; {'time [s]'}]; + if any(strcmp(fieldnames(a),'units')); gdat_data.units=a.units; end + try + [aR,e]=rdaD3D_eff(shot,diag_name,'R-A',exp_name_eff,time_interval); + catch + end + try + [aZ,e]=rdaD3D_eff(shot,diag_name,'z-A',exp_name_eff,time_interval); + catch + disp(['problem with getting z-A in ' diag_name]) + end + if match_rz_to_time + % interpolate R structure on ece data time array, to ease plot vs R + for i=1:length(channels) + radius.data(i,:) = interp1(aR.t,aR.data(channels(i),:),gdat_data.t); + zheight.data(i,:) = interp1(aZ.t,aZ.data(channels(i),:),gdat_data.t); + end + radiuszheight.t=gdat_data.t; + else + radius.data = aR.data(channels,:); + radiuszheight.t=aR.t; + zheight.data = aZ.data(channels,:); + end + ij=find(gdat_data.data<=0); + gdat_data.data(ij)=NaN; + gdat_data.rz.r=radius.data; + gdat_data.rz.z=zheight.data; + gdat_data.rz.t = radiuszheight.t; + gdat_data.data_fullpath = [exp_name_eff '/' diag_name '/Trad-A with R, Z in .r,.z']; + + if strcmp(data_request_eff,'ece_rho') || strcmp(data_request_eff,'eced_rho') + params_equil = gdat_data.gdat_params; + params_equil.data_request = 'equil'; + [equil,params_equil,error_status] = gdat_d3d(shot,params_equil); + if error_status>0 + if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end + return + end + gdat_data.gdat_params.equil = params_equil.equil; + gdat_data.equil = equil; + gdat_data.data_fullpath = [exp_name_eff '/' diag_name '/Trad-A with .r,.z projected on equil ' gdat_data.gdat_params.equil ' in .rhos']; + inb_chord_ece=size(gdat_data.rz.r,1); + inb_time_ece=size(gdat_data.rz.r,2); + psi_out = NaN*ones(inb_chord_ece,inb_time_ece); + rhopolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece); + rhotornorm_out = NaN*ones(inb_chord_ece,inb_time_ece); + rhovolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece); + % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)] + time_equil=[min(gdat_data.rz.t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.rz.t(end)+0.1)]; + for itequil=1:length(time_equil)-1 + rr=equil.Rmesh(:,itequil); + zz=equil.Zmesh(:,itequil); + psirz_in = equil.psi2D(:,:,itequil); + it_ece_inequil = find(gdat_data.rz.t>time_equil(itequil) & gdat_data.rz.t<=time_equil(itequil+1)); + if ~isempty(it_ece_inequil) + r_it_ece_inequil = gdat_data.rz.r(:,it_ece_inequil); + iok = find(~isnan(r_it_ece_inequil) & r_it_ece_inequil>0); + if ~isempty(iok) + rout=r_it_ece_inequil(iok); + z_it_ece_inequil = gdat_data.rz.z(:,it_ece_inequil); + zout=z_it_ece_inequil(iok); + psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout); + [ieff,jeff]=ind2sub(size(gdat_data.rz.r(:,it_ece_inequil)),iok); + for ij=1:length(iok) + psi_out(ieff(ij),it_ece_inequil(jeff(ij))) = psi_at_routzout(ij); + end + rhopolnorm_out(:,it_ece_inequil) = sqrt(abs((psi_out(:,it_ece_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)))); + for it_cx=1:length(it_ece_inequil) + rhotornorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]); + rhovolnorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]); + end + end + end + end + gdat_data.rhos.psi = psi_out; + gdat_data.rhos.rhopolnorm = rhopolnorm_out; + gdat_data.rhos.rhotornorm = rhotornorm_out; + gdat_data.rhos.rhovolnorm = rhovolnorm_out; + gdat_data.rhos.t = gdat_data.rz.t; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'eqdsk'} + % + time_eqdsks=1.; % default time + if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time) + time_eqdsks = gdat_data.gdat_params.time; + else + gdat_data.gdat_params.time = time_eqdsks; + disp(['"time" is expected as an option, choose default time = ' num2str(time_eqdsks)]); + end + gdat_data.gdat_params.time = time_eqdsks; + gdat_data.t = time_eqdsks; + zshift = 0.; + if isfield(gdat_data.gdat_params,'zshift') && ~isempty(gdat_data.gdat_params.zshift) + zshift = gdat_data.gdat_params.zshift; + else + gdat_data.gdat_params.zshift = zshift; + end + gdat_data.gdat_params.zshift = zshift; + params_equil = gdat_data.gdat_params; + params_equil.data_request = 'equil'; + [equil,params_equil,error_status] = gdat_d3d(shot,params_equil); + if error_status>0 + if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end + return + end + gdat_data.gdat_params = params_equil; + gdat_data.equil = equil; + for itime=1:length(time_eqdsks) + time_eff = time_eqdsks(itime); + % use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2 + [eqdskD3D, equil_all_t, equil_t_index]=geteqdskD3D(equil,time_eff,gdat_data.gdat_params.zshift,'source',gdat_data.gdat_params.equil); + eqdskD3D.fnamefull = fullfile(['/tmp/' getenv('USER')],['EQDSK_' num2str(shot) 't' num2str(time_eff,'%.4f')]); + cocos_out = equil.cocos; + if isfield(gdat_data.gdat_params,'cocos') && ~isempty(gdat_data.gdat_params.cocos) + cocos_out = gdat_data.gdat_params.cocos; + else + gdat_data.gdat_params.cocos = cocos_out; + end + if equil.cocos ~= cocos_out + [eqdsk_cocosout, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdskD3D,[equil.cocos cocos_out]); + else + eqdsk_cocosout = eqdskD3D; + end + % for several times, use array of structure for eqdsks, + % cannot use it for psi(R,Z) in .data and .x since R, Z might be different at different times, + % so project psi(R,Z) on Rmesh, Zmesh of 1st time + if length(time_eqdsks) > 1 + gdat_data.eqdsk{itime} = write_eqdsk(eqdskD3D.fnamefull,eqdsk_cocosout); + if itime==1 + gdat_data.data(:,:,itime) = gdat_data.eqdsk{itime}.psi; + gdat_data.dim{1} = gdat_data.eqdsk{itime}.rmesh; + gdat_data.dim{2} = gdat_data.eqdsk{itime}.zmesh; + else + xx=repmat(reshape(gdat_data.dim{1},length(gdat_data.dim{1}),1),1,size(gdat_data.eqdsk{itime}.psi,2)); + yy=repmat(reshape(gdat_data.dim{2},1,length(gdat_data.dim{2})),size(gdat_data.eqdsk{itime}.psi,1),1); + aa = interpos2Dcartesian(gdat_data.eqdsk{itime}.rmesh,gdat_data.eqdsk{itime}.zmesh ... + ,gdat_data.eqdsk{itime}.psi,xx,yy,-1,-1); + gdat_data.data(:,:,itime) = aa; + end + else + gdat_data.eqdsk = write_eqdsk(eqdskD3D.fnamefull,eqdsk_cocosout); + gdat_data.data = gdat_data.eqdsk.psi; + gdat_data.dim{1} = gdat_data.eqdsk.rmesh; + gdat_data.dim{2} = gdat_data.eqdsk.zmesh; + end + end + gdat_data.dim{3} = gdat_data.t; + gdat_data.x = gdat_data.dim(1:2); + gdat_data.data_fullpath=['psi(R,Z) and eqdsk from geteqdskD3D with ' gdat_data.gdat_params.equil ' ;zshift=' num2str(zshift)]; + gdat_data.units = 'T m^2'; + gdat_data.dimunits = {'m','m','s'}; + gdat_data.request_description = ['data=psi, x=(R,Z), eqdsk contains eqdsk structure with which ' ... + 'plot_eqdsk, write_eqdsk, read_eqdsk can be used']; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'equil'} + % get equil params and time array in gdat_data.t + [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data); + % since Lpf depends on time, need to load all first and then loop over time for easier mapping + [qpsi,e]=rdaD3D_eff(shot,DIAG,'Qpsi',exp_name_eff); + ndimrho = size(qpsi.data,2); + if ndimrho==NTIME_Lpf + % data seems to be transposed + ndimrho = size(qpsi.data,1); + itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t + else + itotransposeback = 0; + end + qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback); + ijnan=find(isnan(qpsi.value)); + qpsi.value(ijnan)=0; + [psi_tree,e]=rdaD3D_eff(shot,DIAG,'PFL',exp_name_eff); + psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback); + [phi_tree,e]=rdaD3D_eff(shot,DIAG,'TFLx',exp_name_eff); + phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback); + [Vol,e]=rdaD3D_eff(shot,DIAG,'Vol',exp_name_eff); + Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback); + [Area,e]=rdaD3D_eff(shot,DIAG,'Area',exp_name_eff); + Area=adapt_rda(Area,NTIME,2*ndimrho,itotransposeback); + [Ri,e]=rdaD3D_eff(shot,DIAG,'Ri',exp_name_eff); + Ri=adapt_rda(Ri,NTIME,M_Rmesh,itotransposeback); + [Zj,e]=rdaD3D_eff(shot,DIAG,'Zj',exp_name_eff); + Zj=adapt_rda(Zj,NTIME,N_Zmesh,itotransposeback); + [PFM_tree,e]=rdaD3D_eff(shot,DIAG,'PFM','-raw','-exp',exp_name_eff); % -raw necessary for IDE + PFM_tree=adaptPFM_rda(PFM_tree,M_Rmesh,N_Zmesh,NTIME); + [Pres,e]=rdaD3D_eff(shot,DIAG,'Pres',exp_name_eff); + Pres=adapt_rda(Pres,NTIME,2*ndimrho,itotransposeback); + [Jpol,e]=rdaD3D_eff(shot,DIAG,'Jpol',exp_name_eff); + Jpol=adapt_rda(Jpol,NTIME,2*ndimrho,itotransposeback); + [FFP,e]=rdaD3D_eff(shot,DIAG,'FFP',exp_name_eff); + if ~isempty(FFP.value) + FFP=adapt_rda(FFP,NTIME,ndimrho,itotransposeback); + else + FFP.value=NaN*ones(NTIME,max(Lpf1_t)); + end + if strcmp(DIAG,'EQI') || strcmp(DIAG,'EQH') || strcmp(DIAG,'IDE') + [Rinv,e]=rdaD3D_eff(shot,DIAG,'Rinv',exp_name_eff); + Rinv=adapt_rda(Rinv,NTIME,ndimrho,itotransposeback); + [R2inv,e]=rdaD3D_eff(shot,DIAG,'R2inv',exp_name_eff); + R2inv=adapt_rda(R2inv,NTIME,ndimrho,itotransposeback); + [Bave,e]=rdaD3D_eff(shot,DIAG,'Bave',exp_name_eff); + Bave=adapt_rda(Bave,NTIME,ndimrho,itotransposeback); + [B2ave,e]=rdaD3D_eff(shot,DIAG,'B2ave',exp_name_eff); + B2ave=adapt_rda(B2ave,NTIME,ndimrho,itotransposeback); + if strcmp(DIAG,'IDE') + FTRA.value=[]; + else + [FTRA,e]=rdaD3D_eff(shot,DIAG,'FTRA',exp_name_eff); + FTRA=adapt_rda(FTRA,NTIME,ndimrho,itotransposeback); + end + else + Rinv.value=[]; R2inv.value=[]; Bave.value=[]; B2ave.value=[]; FTRA.value=[]; + end + [LPFx,e]=rdaD3D_eff(shot,DIAG,'LPFx',exp_name_eff); + LPFx.value=LPFx.value(1:NTIME); LPFx.data=LPFx.value; LPFx.t=LPFx.t(1:NTIME); + [PFxx,e]=rdaD3D_eff(shot,DIAG,'PFxx',exp_name_eff); + PFxx=adapt_rda(PFxx,NTIME,max(LPFx.value)+1,itotransposeback); + [RPFx,e]=rdaD3D_eff(shot,DIAG,'RPFx',exp_name_eff); + RPFx=adapt_rda(RPFx,NTIME,max(LPFx.value)+1,itotransposeback); + [zPFx,e]=rdaD3D_eff(shot,DIAG,'zPFx',exp_name_eff); + zPFx=adapt_rda(zPFx,NTIME,max(LPFx.value)+1,itotransposeback); + % seems "LCFS" q-value is far too large, limit to some max (when diverted) + max_qValue = 30.; % Note could just put a NaN on LCFS value since ill-defined when diverted + for it=1:NTIME + Lpf1 = Lpf1_t(it); + % Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part + % change it to (radial,time) and use only Lpf+1 points up to LCFS + ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially + % set NaNs to zeroes + if qpsi.value(ijok(1),1)<0 + gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue); + else + gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue); + end + % get x values + gdat_data.psi(:,it)=psi_tree.value(it,Lpf1:-1:1)'; + gdat_data.psi_axis(it)= gdat_data.psi(1,it); + gdat_data.psi_lcfs(it)= gdat_data.psi(end,it); + gdat_data.rhopolnorm(:,it) = sqrt(abs((gdat_data.psi(:,it)-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it)))); + if strcmp(DIAG,'EQR'); + % q value has only a few values and from center to edge, assume they are from central rhopol values on + % But they are every other point starting from 3rd + ijk=find(gdat_data.qvalue(:,it)~=0); + if length(ijk)>2 + % now shots have non-zero axis values in eqr + rhoeff=gdat_data.rhopolnorm(ijk,it); + qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above + if ijk(1)>1 + rhoeff = [0.; rhoeff]; + qeff = [qeff(1) ;qeff]; + end + ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it))); + qfit = zeros(size(gdat_data.rhopolnorm(:,it))); + qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]); + else + qfit = zeros(size(gdat_data.rhopolnorm(:,it))); + end + gdat_data.qvalue(:,it) = qfit; + end + % get rhotor values + gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)'; + gdat_data.rhotornorm(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.phi(end,it))); + % get rhovol values + gdat_data.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)'; + gdat_data.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi + gdat_data.rhovolnorm(:,it) = sqrt(abs(gdat_data.vol(:,it) ./ gdat_data.vol(end,it))); + gdat_data.area(:,it)=Area.value(it,2*Lpf1-1:-2:1)'; + gdat_data.dareadpsi(:,it)=Area.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi + gdat_data.Rmesh(:,it) = Ri.value(it,1:M_Rmesh); + gdat_data.Zmesh(:,it) = Zj.value(it,1:N_Zmesh); + gdat_data.psi2D(1:M_Rmesh,1:N_Zmesh,it) = PFM_tree.value(1:M_Rmesh,1:N_Zmesh,it); + gdat_data.pressure(:,it)=Pres.value(it,2*Lpf1-1:-2:1)'; + gdat_data.dpressuredpsi(:,it)=Pres.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi + if ~isempty(Jpol.value) + gdat_data.jpol(:,it)=Jpol.value(it,2*Lpf1-1:-2:1)'; + gdat_data.djpolpsi(:,it)=Jpol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi + else + gdat_data.jpol = []; + gdat_data.djpolpsi = []; + end + gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)'; + gdat_data.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1); + gdat_data.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1); + gdat_data.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1); + if ~isempty(Rinv.value) + gdat_data.rinv(:,it) = Rinv.value(it,Lpf1:-1:1)'; + else + gdat_data.rinv = []; + end + if ~isempty(R2inv.value) + gdat_data.r2inv(:,it) = R2inv.value(it,Lpf1:-1:1)'; + else + gdat_data.r2inv = []; + end + if ~isempty(Bave.value) + gdat_data.bave(:,it) = Bave.value(it,Lpf1:-1:1)'; + else + gdat_data.bave = []; + end + if ~isempty(B2ave.value) + gdat_data.b2ave(:,it) = B2ave.value(it,Lpf1:-1:1)'; + else + gdat_data.b2ave = []; + end + if ~isempty(FTRA.value) + gdat_data.ftra(:,it) = FTRA.value(it,Lpf1:-1:1)'; + else + gdat_data.ftra = []; + end + % + end + gdat_data.x = gdat_data.rhopolnorm; + % + [equil_Rcoil,e]=rdaD3D_eff(shot,DIAG,'Rcl',exp_name_eff); + gdat_data.Rcoils=equil_Rcoil.value; + [equil_Zcoil,e]=rdaD3D_eff(shot,DIAG,'Zcl',exp_name_eff); + gdat_data.Zcoils=equil_Zcoil.value; + % + if strcmp(DIAG,'IDE') + [IpiPSI,e]=rdaD3D_eff(shot,'IDG','Itor',exp_name_eff); + if length(IpiPSI.value)~=NTIME + disp('Problem with nb time points of IDG/Itor with respect to IDE NTIME, check with O. Sauter') + return + end + else + [IpiPSI,e]=rdaD3D_eff(shot,DIAG,'IpiPSI',exp_name_eff); + end + gdat_data.Ip = IpiPSI.value(1:NTIME); + % + gdat_data.data = gdat_data.qvalue; % put q in data + gdat_data.units=[]; % not applicable + gdat_data.data_fullpath = [DIAG ' from expname: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)']; + gdat_data.cocos = 17; % should check FPP + gdat_data.dim{1} = gdat_data.x; + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits{1} = 'rhopolnorm'; + gdat_data.dimunits{2} = 'time [s]'; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'ne','te'} + exp_name_eff = gdat_data.gdat_params.exp_name; + % ne or Te from Thomson data on raw z mesh vs (z,t) + nodenameeff = [upper(data_request_eff(1)) 'e_c']; + node_child_nameeff = [upper(data_request_eff(1)) 'e_core']; + [a,error_status]=rdaD3D_eff(shot,'VTA',nodenameeff,exp_name_eff); + if isempty(a.data) || isempty(a.t) || error_status>0 + if gdat_params.nverbose>=3; + a + disp(['with data_request= ' data_request_eff]) + end + return + end + gdat_data.(lower(node_child_nameeff)).data = a.data; + gdat_data.(lower(node_child_nameeff)).t = a.t; + gdat_data.t = a.t; + if any(strcmp(fieldnames(a),'units')) + gdat_data.units=a.units; + end + [alow,e]=rdaD3D_eff(shot,'VTA',[upper(data_request_eff(1)) 'elow_c'],exp_name_eff); + [aup,e]=rdaD3D_eff(shot,'VTA',[upper(data_request_eff(1)) 'eupp_c'],exp_name_eff); + gdat_data.(lower(node_child_nameeff)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff)).data,gdat_data.(lower(node_child_nameeff)).data-alow.value); + [a,error_status]=rdaD3D_eff(shot,'VTA','R_core',exp_name_eff); + gdat_data.(lower(node_child_nameeff)).r = repmat(a.data,size(gdat_data.(lower(node_child_nameeff)).data,1),1); + [a,error_status]=rdaD3D_eff(shot,'VTA','Z_core',exp_name_eff); + gdat_data.(lower(node_child_nameeff)).z = repmat(a.data',1,size(gdat_data.(lower(node_child_nameeff)).data,2)); + gdat_data.dimunits=[{'Z [m]'} ; {'time [s]'}]; + gdat_data.data_fullpath=[data_request_eff ' from VTA/' upper(data_request_eff(1)) 'e_c and ' upper(data_request_eff(1)) 'e_e']; + nb_core = size(gdat_data.(lower(node_child_nameeff)).z,1); + gdat_data.data = []; + gdat_data.data(1:nb_core,:) = gdat_data.(lower(node_child_nameeff)).data(1:nb_core,:); + gdat_data.dim=[{gdat_data.(lower(node_child_nameeff)).z};{gdat_data.(lower(node_child_nameeff)).t}]; + gdat_data.error_bar(1:nb_core,:) = gdat_data.(lower(node_child_nameeff)).error_bar(1:nb_core,:); + + % add edge part + nodenameeff_e = [upper(data_request_eff(1)) 'e_e']; + node_child_nameeff_e = [upper(data_request_eff(1)) 'e_edge']; + [a,error_status]=rdaD3D_eff(shot,'VTA',nodenameeff_e,exp_name_eff); + gdat_data.(lower(node_child_nameeff_e)).data = a.data; + ij_edge_notok = find(a.data>5e21); + gdat_data.(lower(node_child_nameeff_e)).data(ij_edge_notok) = NaN; + gdat_data.(lower(node_child_nameeff_e)).t = a.t; + if ~isempty(a.data) + [a,error_status]=rdaD3D_eff(shot,'VTA','R_edge',exp_name_eff); + gdat_data.(lower(node_child_nameeff_e)).r = repmat(a.data,size(gdat_data.(lower(node_child_nameeff_e)).data,1),1); + [a,error_status]=rdaD3D_eff(shot,'VTA','Z_edge',exp_name_eff); + gdat_data.(lower(node_child_nameeff_e)).z = repmat(a.data',1,size(gdat_data.(lower(node_child_nameeff_e)).data,2)); + nb_edge = size(gdat_data.(lower(node_child_nameeff_e)).z,1); + iaaa=iround_os(gdat_data.(lower(node_child_nameeff_e)).t,gdat_data.(lower(node_child_nameeff)).t); + gdat_data.data(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).data(1:nb_edge,iaaa); + gdat_data.dim{1}(nb_core+1:nb_core+nb_edge,:)=gdat_data.(lower(node_child_nameeff_e)).z(1:nb_edge,iaaa); + [alow,e]=rdaD3D_eff(shot,'VTA',[upper(data_request_eff(1)) 'elow_e'],exp_name_eff); + [aup,e]=rdaD3D_eff(shot,'VTA',[upper(data_request_eff(1)) 'eupp_e'],exp_name_eff); + gdat_data.(lower(node_child_nameeff_e)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff_e)).data, ... + gdat_data.(lower(node_child_nameeff_e)).data-alow.value); + gdat_data.error_bar(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).error_bar(1:nb_edge,iaaa); + else + gdat_data.(lower(node_child_nameeff_e)).data = []; + gdat_data.(lower(node_child_nameeff_e)).t = []; + gdat_data.(lower(node_child_nameeff_e)).error_bar = []; + gdat_data.(lower(node_child_nameeff_e)).r = []; + gdat_data.(lower(node_child_nameeff_e)).z = []; + end + gdat_data.x=gdat_data.dim{1}; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'ne_rho', 'te_rho', 'nete_rho'} + if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit) || ~isnumeric(gdat_data.gdat_params.fit) + gdat_data.gdat_params.fit = 1; % default do fit + end + params_eff = gdat_data.gdat_params; + params_eff.data_request=data_request_eff(1:2); % start with ne if nete_rho + % get raw data + [gdat_data,params_kin,error_status]=gdat_d3d(shot,params_eff); + gdat_data.gdat_params.data_request=data_request_eff; + gdat_data.gdat_request=data_request_eff; + if error_status>0 + if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end + return + end + % add rho coordinates + params_eff.data_request='equil'; + [equil,params_equil,error_status]=gdat_d3d(shot,params_eff); + if error_status>0 + if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end + return + end + gdat_data.gdat_params.equil = params_equil.equil; + %gdat_data.equil = equil; + + % core + node_child_nameeff = [upper(data_request_eff(1)) 'e_core']; + inb_chord_core=size(gdat_data.(lower(node_child_nameeff)).r,1); + inb_time_core=size(gdat_data.(lower(node_child_nameeff)).r,2); + psi_out_core = NaN*ones(inb_chord_core,inb_time_core); + rhopolnorm_out_core = NaN*ones(inb_chord_core,inb_time_core); + rhotornorm_out_core = NaN*ones(inb_chord_core,inb_time_core); + rhovolnorm_out_core = NaN*ones(inb_chord_core,inb_time_core); + % edge + node_child_nameeff_e = [upper(data_request_eff(1)) 'e_edge']; + inb_chord_edge=size(gdat_data.(lower(node_child_nameeff_e)).r,1); + inb_time_edge=size(gdat_data.(lower(node_child_nameeff_e)).r,2); + psi_out_edge = NaN*ones(inb_chord_edge,inb_time_edge); + rhopolnorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge); + rhotornorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge); + rhovolnorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge); + % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)] + time_equil=[min(gdat_data.(lower(node_child_nameeff)).t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.(lower(node_child_nameeff)).t(end)+0.1)]; + for itequil=1:length(time_equil)-1 + rr=equil.Rmesh(:,itequil); + zz=equil.Zmesh(:,itequil); + psirz_in = equil.psi2D(:,:,itequil); + it_core_inequil = find(gdat_data.(lower(node_child_nameeff)).t>time_equil(itequil) & gdat_data.(lower(node_child_nameeff)).t<=time_equil(itequil+1)); + if ~isempty(it_core_inequil) + rout_core=gdat_data.(lower(node_child_nameeff)).r(:,it_core_inequil); + zout_core=gdat_data.(lower(node_child_nameeff)).z(:,it_core_inequil); + psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core); + psi_out_core(:,it_core_inequil) = reshape(psi_at_routzout,inb_chord_core,length(it_core_inequil)); + rhopolnorm_out_core(:,it_core_inequil) = sqrt((psi_out_core(:,it_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))); + for it_cx=1:length(it_core_inequil) + rhotornorm_out_core(:,it_core_inequil(it_cx)) = ... + interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]); + rhovolnorm_out_core(:,it_core_inequil(it_cx)) = ... + interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]); + end + end + % edge + it_edge_inequil = find(gdat_data.(lower(node_child_nameeff_e)).t>time_equil(itequil) & gdat_data.(lower(node_child_nameeff_e)).t<=time_equil(itequil+1)); + if ~isempty(it_edge_inequil) + rout_edge=gdat_data.(lower(node_child_nameeff_e)).r(:,it_edge_inequil); + zout_edge=gdat_data.(lower(node_child_nameeff_e)).z(:,it_edge_inequil); + psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_edge,zout_edge); + psi_out_edge(:,it_edge_inequil) = reshape(psi_at_routzout,inb_chord_edge,length(it_edge_inequil)); + rhopolnorm_out_edge(:,it_edge_inequil) = sqrt((psi_out_edge(:,it_edge_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))); + for it_cx=1:length(it_edge_inequil) + rhotornorm_out_edge(:,it_edge_inequil(it_cx)) = ... + interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]); + rhovolnorm_out_edge(:,it_edge_inequil(it_cx)) = ... + interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]); + end + end + end + gdat_data.(lower(node_child_nameeff)).psi = psi_out_core; + gdat_data.(lower(node_child_nameeff)).rhopolnorm = rhopolnorm_out_core; + gdat_data.(lower(node_child_nameeff)).rhotornorm = rhotornorm_out_core; + gdat_data.(lower(node_child_nameeff)).rhovolnorm = rhovolnorm_out_core; + gdat_data.(lower(node_child_nameeff_e)).psi = psi_out_edge; + gdat_data.(lower(node_child_nameeff_e)).rhopolnorm = rhopolnorm_out_edge; + gdat_data.(lower(node_child_nameeff_e)).rhotornorm = rhotornorm_out_edge; + gdat_data.(lower(node_child_nameeff_e)).rhovolnorm = rhovolnorm_out_edge; + % put values of rhopolnorm for dim{1} by default, all radial mesh for combined core, edge in grids_1d + gdat_data.x = gdat_data.(lower(node_child_nameeff)).rhopolnorm; + iaaa=iround_os(gdat_data.(lower(node_child_nameeff_e)).t,gdat_data.(lower(node_child_nameeff)).t); + gdat_data.x(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhopolnorm(:,iaaa); + gdat_data.dim{1} = gdat_data.x; + gdat_data.dimunits{1} = 'rhopolnorm'; + gdat_data.grids_1d.rhopolnorm = gdat_data.x; + gdat_data.grids_1d.psi = gdat_data.(lower(node_child_nameeff)).psi; + gdat_data.grids_1d.psi(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).psi(:,iaaa); + gdat_data.grids_1d.rhotornorm = gdat_data.(lower(node_child_nameeff)).rhotornorm; + gdat_data.grids_1d.rhotornorm(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhotornorm(:,iaaa); + gdat_data.grids_1d.rhovolnorm = gdat_data.(lower(node_child_nameeff)).rhovolnorm; + gdat_data.grids_1d.rhovolnorm(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhovolnorm(:,iaaa); + + gdat_data.data_fullpath = [gdat_data.data_fullpath ' projected on equilibrium ' gdat_data.gdat_params.equil]; + + % if nete_rho, copy data as .ne, get .te and put pe=e ne Te in data: + gdat_data.(data_request_eff(1:2)).data = gdat_data.data; + gdat_data.(data_request_eff(1:2)).error_bar = gdat_data.error_bar; + gdat_data.(data_request_eff(1:2)).units = gdat_data.units; + gdat_data.(data_request_eff(1:2)).core = gdat_data.([data_request_eff(1:2) '_core']); + gdat_data.(data_request_eff(1:2)).edge = gdat_data.([data_request_eff(1:2) '_edge']); + gdat_data = rmfield(gdat_data,{[data_request_eff(1:2) '_core'],[data_request_eff(1:2) '_edge']}); + if strcmp(data_request_eff(1:4),'nete') + params_eff.data_request=data_request_eff(3:4); + [gdat_data_te,params_kin,error_status]=gdat_d3d(shot,params_eff); + gdat_data.te.data = gdat_data_te.data; + gdat_data.te.error_bar = gdat_data_te.error_bar; + gdat_data.te.units = gdat_data_te.units; + gdat_data.te.core = gdat_data_te.te_core; + gdat_data.te.edge = gdat_data_te.te_edge; + gdat_data.te.error_bar = gdat_data_te.error_bar; + gdat_data.te.core.psi = gdat_data.ne.core.psi; + gdat_data.te.core.rhopolnorm = gdat_data.ne.core.rhopolnorm; + gdat_data.te.core.rhotornorm = gdat_data.ne.core.rhotornorm; + gdat_data.te.core.rhovolnorm = gdat_data.ne.core.rhovolnorm; + gdat_data.te.edge.psi = gdat_data.ne.edge.psi; + gdat_data.te.edge.rhopolnorm = gdat_data.ne.edge.rhopolnorm; + gdat_data.te.edge.rhotornorm = gdat_data.ne.edge.rhotornorm; + gdat_data.te.edge.rhovolnorm = gdat_data.ne.edge.rhovolnorm; + % put pe in main gdat_data + gdat_data.data = 1.6022e-19.*gdat_data.ne.data.*gdat_data.te.data; + gdat_data.error_bar = 1.6022e-19 .* (gdat_data.ne.data .* gdat_data.te.error_bar ... + + gdat_data.te.data .* gdat_data.ne.error_bar); + gdat_data.units='N/m^2; 1.6022e-19 ne Te'; + gdat_data.data_fullpath=['pe=1.6e-19*ne*Te in data, .ne, .te ' gdat_data.data_fullpath(3:end)]; + gdat_data.label = 'pe'; + end + + % defaults for fits, so user always gets std structure + gdat_data.fit.rhotornorm = []; % same for both ne and te + gdat_data.fit.rhopolnorm = []; + gdat_data.fit.t = []; + gdat_data.fit.te.data = []; + gdat_data.fit.te.drhotornorm = []; + gdat_data.fit.ne.data = []; + gdat_data.fit.ne.drhotornorm = []; + gdat_data.fit.raw.te.data = []; + gdat_data.fit.raw.te.rhotornorm = []; + gdat_data.fit.raw.ne.data = []; + gdat_data.fit.raw.ne.rhotornorm = []; + fit_tension_default = -0.1; + if isfield(gdat_data.gdat_params,'fit_tension') + fit_tension = gdat_data.gdat_params.fit_tension; + else + fit_tension = fit_tension_default; + end + if ~isstruct(fit_tension) + fit_tension_eff.te = fit_tension; + fit_tension_eff.ne = fit_tension; + fit_tension = fit_tension_eff; + else + if ~isfield(fit_tension,'te'); fit_tension.te = fit_tension_default; end + if ~isfield(fit_tension,'ne '); fit_tension.ne = fit_tension_default; end + end + gdat_data.gdat_params.fit_tension = fit_tension; + if isfield(gdat_data.gdat_params,'fit_nb_rho_points') + fit_nb_rho_points = gdat_data.gdat_params.fit_nb_rho_points; + else + fit_nb_rho_points = 201; + end + gdat_data.gdat_params.fit_nb_rho_points = fit_nb_rho_points; + % + if gdat_data.gdat_params.fit==1 + % add fits + gdat_data.fit.t = gdat_data.t; + rhotornormfit = linspace(0,1,fit_nb_rho_points)'; + gdat_data.fit.rhotornorm = rhotornormfit; + gdat_data.fit.rhopolnorm = NaN*ones(length(rhotornormfit),length(gdat_data.fit.t)); + if any(strfind(data_request_eff,'te')) + gdat_data.fit.raw.te.data = NaN*ones(size(gdat_data.te.data)); + gdat_data.fit.raw.te.error_bar = NaN*ones(size(gdat_data.te.data)); + gdat_data.fit.raw.te.rhotornorm = NaN*ones(size(gdat_data.te.data)); + gdat_data.fit.te.data = gdat_data.fit.rhopolnorm; + gdat_data.fit.te.drhotornorm = gdat_data.fit.rhopolnorm; + end + if any(strfind(data_request_eff,'ne')) + gdat_data.fit.raw.ne.data = NaN*ones(size(gdat_data.ne.data)); + gdat_data.fit.raw.ne.error_bar = NaN*ones(size(gdat_data.ne.data)); + gdat_data.fit.raw.ne.rhotornorm = NaN*ones(size(gdat_data.ne.data)); + gdat_data.fit.ne.data =gdat_data.fit.rhopolnorm; + gdat_data.fit.ne.drhotornorm = gdat_data.fit.rhopolnorm; + end + for it=1:length(gdat_data.t) + % make rhotor->rhopol transformation for each time since equilibrium might have changed + irhook=find(gdat_data.grids_1d.rhotornorm(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<1); % no need for ~isnan + [rhoeff isort]=sort(gdat_data.grids_1d.rhotornorm(irhook,it)); + gdat_data.fit.rhopolnorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.grids_1d.rhopolnorm(irhook(isort),it); 1],rhotornormfit,-0.1,[2 2],[0 1]); + if any(strfind(data_request_eff,'te')) + idatate = find(gdat_data.te.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05); + if length(idatate)>0 + gdat_data.fit.raw.te.rhotornorm(idatate,it) = gdat_data.grids_1d.rhotornorm(idatate,it); + gdat_data.fit.raw.te.data(idatate,it) = gdat_data.te.data(idatate,it); + gdat_data.fit.raw.te.error_bar(idatate,it) = gdat_data.te.error_bar(idatate,it); + [rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhotornorm(idatate,it)); + rhoeff = [0; rhoeff]; + teeff = gdat_data.te.data(idatate(irhoeff),it); + te_err_eff = gdat_data.te.error_bar(idatate(irhoeff),it); + % they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar + ij=find(teeff./te_err_eff>10.); + if ~isempty(ij); te_err_eff(ij) = teeff(ij)./0.1; end + % + teeff = [teeff(1); teeff]; + te_err_eff = [1e4; te_err_eff]; + [gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhotornorm(:,it)] = interpos(rhoeff,teeff,rhotornormfit,fit_tension.te,[1 0],[0 0],te_err_eff); + end + end + if any(strfind(data_request_eff,'ne')) + idatane = find(gdat_data.ne.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05); + if length(idatane)>0 + gdat_data.fit.raw.ne.rhotornorm(idatane,it) = gdat_data.grids_1d.rhotornorm(idatane,it); + gdat_data.fit.raw.ne.data(idatane,it) = gdat_data.ne.data(idatane,it); + gdat_data.fit.raw.ne.error_bar(idatane,it) = gdat_data.ne.error_bar(idatane,it); + [rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhotornorm(idatane,it)); + rhoeff = [0; rhoeff]; + % they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar + neeff = gdat_data.ne.data(idatane(irhoeff),it); + ne_err_eff = gdat_data.ne.error_bar(idatane(irhoeff),it); + ij=find(neeff./ne_err_eff>10.); + if ~isempty(ij); ne_err_eff(ij) = neeff(ij)./0.1; end + % + neeff = [neeff(1); neeff]; + ne_err_eff = [1e21; ne_err_eff]; + [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhotornorm(:,it)] = interpos(rhoeff,neeff,rhotornormfit,fit_tension.ne,[1 0],[0 0],ne_err_eff); + end + end + end + end + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'pgyro'} + % LOAD MULTI CHANNEL DATA ECS + % powers, frequencies, etc + params_eff = gdat_data.gdat_params; + params_eff.data_request={'ECS','PECRH'}; + % pgyro tot in index=9 + try + gdat_data=gdat_d3d(shot,params_eff); + gdat_data.data_request = data_request_eff; + gdat_data.gdat_params.data_request = data_request_eff; + catch + if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end + gdat_data.data_request = data_request_eff; + return + end + nb_timepoints = length(gdat_data.t); + pgyro = NaN*ones(nb_timepoints,9); + pgyro(:,9) = reshape(gdat_data.data,nb_timepoints,1); + gdat_data.data = pgyro; + labels{9} = 'ECtot'; + for i=1:4 + % "old" ECRH1 gyrotrons: gyro 1 to 4 in pgyro + params_eff.data_request={'ECS',['PG' num2str(i)]}; + gdat_data_i=gdat_d3d(shot,params_eff); + if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) + else + gdat_data.ec{i} = gdat_data_i; + gdat_data.data(:,i) = reshape(gdat_data_i.data,nb_timepoints,1); + gdat_data.dim = [{gdat_data_i.t} {[1:9]}]; + if max(gdat_data_i.data) > 0. + labels{i} = ['EC_' num2str(i)]; + end + end + try + a = sf2par('ECS',shot,'gyr_freq',['P_sy1_g' num2str(i)]); + catch + % gyr_freq not present (old shots for example) + a=[]; + end + if isempty(a) + else + gdat_data.ec{i}.freq = a; + gdat_data.freq_ec(i) = a.value; + end + try + a = sf2par('ECS',shot,'GPolPos',['P_sy1_g' num2str(i)]); + catch + % GPolPos not present + a=[]; + end + if isempty(a) + else + gdat_data.ec{i}.polpos = a.value; + gdat_data.polpos_ec(i) = a.value; + end + try + a = sf2par('ECS',shot,'GTorPos',['P_sy1_g' num2str(i)]); + catch + a=[]; + end + if isempty(a) + else + gdat_data.ec{i}.torpos = a.value; + gdat_data.torpos_ec(i) = a.value; + end + % "new" ECRH2 gyrotrons: gyro 5 to 8 in pgyro + params_eff.data_request={'ECS',['PG' num2str(i) 'N']}; + gdat_data_i=gdat_d3d(shot,params_eff); + if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) + else + gdat_data.ec{i+4} = gdat_data_i; + gdat_data.data(:,i+4) = reshape(gdat_data_i.data,nb_timepoints,1); + gdat_data.dim = [{gdat_data_i.t} {[1:9]}]; + if max(gdat_data_i.data) > 0. + labels{i+4} = ['EC_' num2str(i+4)]; + end + end + try + a = sf2par('ECS',shot,'gyr_freq',['P_sy2_g' num2str(i)]); + catch + a=[]; + end + if isempty(a) + else + gdat_data.ec{i+4}.freq = a; + gdat_data.freq_ec(i+4) = a.value; + end + try + a = sf2par('ECS',shot,'GPolPos',['P_sy2_g' num2str(i)]); + catch + a=[]; + end + if isempty(a) + else + gdat_data.ec{i+4}.polpos = a.value; + gdat_data.polpos_ec(i+4) = a.value; + end + try + a = sf2par('ECS',shot,'GTorPos',['P_sy2_g' num2str(i)]); + catch + a=[]; + end + if isempty(a) + else + gdat_data.ec{i+4}.torpos = a.value; + gdat_data.torpos_ec(i+4) = a.value; + end + params_eff.data_request={'ECN',['G' num2str(i) 'POL']}; + gdat_data_i=gdat_d3d(shot,params_eff); + if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) + else + gdat_data.ec{i+4}.gpol_ec = gdat_data_i; + end + params_eff.data_request={'ECN',['G' num2str(i) 'TOR']}; + gdat_data_i=gdat_d3d(shot,params_eff); + if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) + else + gdat_data.ec{i+4}.gtor_ec = gdat_data_i; + end + params_eff.data_request={'ECN',['G' num2str(i) 'PO4']}; + gdat_data_i=gdat_d3d(shot,params_eff); + if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) + else + gdat_data.ec{i+4}.gpo4_ec = gdat_data_i; + end + params_eff.data_request={'ECN',['G' num2str(i) 'PO8']}; + gdat_data_i=gdat_d3d(shot,params_eff); + if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) + else + gdat_data.ec{i+4}.gpo8_ec = gdat_data_i; + end + end + if ~isempty(gdat_data.dim) + gdat_data.t = gdat_data.dim{1}; + gdat_data.x = gdat_data.dim{2}; + gdat_data.dimunits=[{'time [s]'} {'ECRH1(1:4) ECRH2(1:4) ECtot'}]; + gdat_data.units='W'; + gdat_data.freq_ech_units = 'GHz'; + gdat_data.data_fullpath=['ECS/' 'PGi and PGiN, etc']; + icount=0; + for i=1:length(labels) + if ~isempty(labels{i}) + icount=icount+1; + gdat_data.label{icount} = labels{i}; + end + end + else + gdat_data.freq_ech_units =[]'; + gdat_data.ec = []; + gdat_data.freq_ec = []; + gdat_data.polpos_ec = []; + gdat_data.torpos_ec = []; + end + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'powers'} + sources_avail = {'ohm','ec','nbi','ic'}; % note should allow ech, nbh, ohmic in parameter sources + for i=1:length(sources_avail) + gdat_data.(sources_avail{i}).data = []; + gdat_data.(sources_avail{i}).units = []; + gdat_data.(sources_avail{i}).dim=[]; + gdat_data.(sources_avail{i}).dimunits=[]; + gdat_data.(sources_avail{i}).t=[]; + gdat_data.(sources_avail{i}).x=[]; + gdat_data.(sources_avail{i}).data_fullpath=[]; + gdat_data.(sources_avail{i}).label=[]; + end + if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) + gdat_data.gdat_params.source = sources_avail; + elseif ~iscell(gdat_data.gdat_params.source) + if ischar(gdat_data.gdat_params.source) + gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source); + if ~any(strmatch(gdat_data.gdat_params.source,lower(sources_avail))) + if (gdat_params.nverbose>=1) + warning(['source= ' gdat_data.gdat_params.source ' is not part of the available sources: ' sprintf('''%s'' ',sources_avail{:})]); + end + return + else + gdat_data.gdat_params.source = {gdat_data.gdat_params.source}; + end + else + if (gdat_params.nverbose>=1); warning([' source parameter not compatible with: ' sprintf('''%s'' ',sources_avail{:})]); end + return + end + else + for i=1:length(gdat_data.gdat_params.source) + gdat_data.gdat_params.source{i} = lower(gdat_data.gdat_params.source{i}); + if ~any(strmatch(gdat_data.gdat_params.source{i},lower(sources_avail))) + if gdat_data.gdat_params.nverbose>=1 + warning(['source = ' gdat_data.gdat_params.source{i} ' not expected with data_request= ' data_request_eff]) + end + end + end + end + % always start from ohmic so can use this time as base time since should yield full shot + + fields_to_copy = {'data','units','dim','dimunits','t','x','data_fullpath','label','help','gdat_params'}; + fields_to_not_copy = {'shot','gdat_request'}; + % total of each source in .data, but full data in subfield like pgyro in .ec, to check for nbi + params_eff = gdat_data.gdat_params; + % ohmic, use its time-base + params_eff.data_request={'TOT','P_OH'}; + try + ohm=gdat_d3d(shot,params_eff); + catch + ohm.data = []; + ohm.dim = []; + end + if ~isempty(ohm.data) && ~isempty(ohm.dim) + for i=1:length(fields_to_copy) + if isfield(ohm,fields_to_copy{i}) + gdat_data.ohm.(fields_to_copy{i}) = ohm.(fields_to_copy{i}); + end + end + gdat_data.ohm.raw_data = gdat_data.ohm.data; + else + if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end + return + end + mapping_for_d3d.timedim = 1; mapping_for_d3d.gdat_timedim = 1; + taus = -10; + % + % add each source in main.data, on ohm time array + gdat_data.t = gdat_data.ohm.t; + gdat_data.dim{1} = gdat_data.t; + gdat_data.dimunits{1} = 's'; + gdat_data.ohm.data = interpos(gdat_data.t,gdat_data.ohm.raw_data,5.*taus); + gdat_data.data = reshape(gdat_data.ohm.data,length(gdat_data.t),1); + gdat_data.ohm.tension = 5.*taus; + gdat_data.x =[1]; + gdat_data.label={'P_{ohm}'}; + gdat_data.units = 'W'; + % + if any(strmatch('ec',gdat_data.gdat_params.source)) + % ec + params_eff.data_request={'ECS','PECRH'}; + params_eff.data_request='pgyro'; + try + ec=gdat_d3d(shot,params_eff); + catch + end + if ~isempty(ec.data) && ~isempty(ec.dim) + for i=1:length(fields_to_copy) + % if has pgyro, use not_copy + if isfield(ec,fields_to_copy{i}) && ~any(strmatch(fields_to_not_copy,fields_to_copy{i})) + gdat_data.ec.(fields_to_copy{i}) = ec.(fields_to_copy{i}); + end + end + gdat_data.data(:,end+1) = interpos(-21,gdat_data.ec.t,gdat_data.ec.data(:,end),gdat_data.t); + gdat_data.x(end+1) =gdat_data.x(end)+1; + gdat_data.label{end+1}='P_{ec}'; + end + end + % + if any(strmatch('nb',gdat_data.gdat_params.source)) + % nbi + params_eff.data_request={'NIS','PNIQ'}; + try + nbi=gdat_d3d(shot,params_eff); + catch + end + if ~isempty(nbi.data) && ~isempty(nbi.dim) + for i=1:length(fields_to_copy) + if isfield(nbi,fields_to_copy{i}) + gdat_data.nbi.(fields_to_copy{i}) = nbi.(fields_to_copy{i}); + end + end + % add to main with linear interpolation and 0 for extrapolated values + gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi.t,gdat_data.nbi.data(:,end),gdat_data.t); + gdat_data.x(end+1) =gdat_data.x(end)+1; + gdat_data.label{end+1}='P_{nbi}'; + end + end + % + if any(strmatch('ic',gdat_data.gdat_params.source)) + % ic + params_eff.data_request={'ICP','PICRN'}; + try + ic=gdat_d3d(shot,params_eff); + catch + end + if ~isempty(ic.data) && ~isempty(ic.dim) + for i=1:length(fields_to_copy) + if isfield(ic,fields_to_copy{i}) + gdat_data.ic.(fields_to_copy{i}) = ic.(fields_to_copy{i}); + end + end + gdat_data.data(:,end+1) = interpos(-21,gdat_data.ic.t,gdat_data.ic.data,gdat_data.t); + gdat_data.x(end+1) =gdat_data.x(end)+1; + gdat_data.label{end+1}='P_{ic}'; + end + end + % add tot power + gdat_data.data(:,end+1) = sum(gdat_data.data,2); + gdat_data.label{end+1}='P_{tot}'; + gdat_data.x(end+1) =gdat_data.x(end)+1; + gdat_data.dim{2} = gdat_data.x; gdat_data.dimunits{2} = ''; + gdat_data.data_fullpath = 'tot powers from each sources, and total power in .data(:,end)'; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'q_rho'} + [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data); + % since Lpf depends on time, need to load all first and then loop over time for easier mapping + [qpsi,e]=rdaD3D_eff(shot,DIAG,'Qpsi',exp_name_eff); + ndimrho = size(qpsi.data,2); + if ndimrho==NTIME_Lpf + % data seems to be transposed + ndimrho = size(qpsi.data,1); + itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t + else + itotransposeback = 0; + end + qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback); + ijnan=find(isnan(qpsi.value)); + qpsi.value(ijnan)=0; + [psi_tree,e]=rdaD3D_eff(shot,DIAG,'PFL',exp_name_eff); + psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback); + [phi_tree,e]=rdaD3D_eff(shot,DIAG,'TFLx',exp_name_eff); + phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback); + [Vol,e]=rdaD3D_eff(shot,DIAG,'Vol',exp_name_eff); + Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback); + % seems "LCFS" q-value is far too large, limit to some max (when diverted) + max_qValue = 30.; % Note could just put a NaN on LCFS value since ill-defined when diverted + for it=1:NTIME + Lpf1 = Lpf1_t(it); + % Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part + % change it to (radial,time) and use only Lpf+1 points up to LCFS + ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially + % set NaNs to zeroes + if qpsi.value(ijok(1),1)<0 + gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue); + else + gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue); + end + % get x values + psi_it=psi_tree.value(it,Lpf1:-1:1)'; + gdat_data.psi_axis(it)= psi_it(1); + gdat_data.psi_lcfs(it)= psi_it(end); + gdat_data.rhopolnorm(:,it) = sqrt(abs((psi_it-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it)))); + if strcmp(DIAG,'EQR'); + % q value has only a few values and from center to edge, assume they are from central rhopol values on + % But they are every other point starting from 3rd + ijk=find(gdat_data.qvalue(:,it)~=0); + if length(ijk)>2 + % now shots have non-zero axis values in eqr + rhoeff=gdat_data.rhopolnorm(ijk,it); + qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above + if ijk(1)>1 + rhoeff = [0.; rhoeff]; + qeff = [qeff(1) ;qeff]; + end + ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it))); + qfit = zeros(size(gdat_data.rhopolnorm(:,it))); + qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]); + else + qfit = zeros(size(gdat_data.rhopolnorm(:,it))); + end + gdat_data.qvalue(:,it) = qfit; + end + % get rhotor values + phi_it = phi_tree.value(it,Lpf1:-1:1)'; + gdat_data.rhotornorm(:,it) = sqrt(abs(phi_it ./ phi_it(end))); + % get rhovol values + vol_it=Vol.value(it,2*Lpf1-1:-2:1)'; + gdat_data.rhovolnorm(:,it) = sqrt(abs(vol_it ./ vol_it(end))); + % Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part + % change it to (radial,time) and use only Lpf+1 points up to LCFS + ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially + % max value 1e6, change to 2*extrapolation (was max_qValue before) + q_edge=interpos(gdat_data.rhotornorm(1:end-1,it),qpsi.value(it,Lpf1:-1:2),1,-0.1); + gdat_data.qvalue(:,it) = qpsi.value(it,Lpf1:-1:1)'; + if abs(gdat_data.qvalue(end,it)) > 1e3 + % assume diverted + gdat_data.qvalue(end,it) = 2. * q_edge; + end +% $$$ if qpsi.value(ijok(1),1)<0 +% $$$ gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue); +% $$$ else +% $$$ gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue); +% $$$ end + end + gdat_data.x = gdat_data.rhopolnorm; + % get time values + gdat_data.data = gdat_data.qvalue; % put q in data + gdat_data.units=[]; % not applicable + gdat_data.data_fullpath = [DIAG ' from expname: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)']; + gdat_data.cocos = 17; % should check FPP + gdat_data.dim{1} = gdat_data.x; + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits{1} = 'rhopolnorm'; + gdat_data.dimunits{2} = 'time [s]'; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'psi_axis', 'psi_edge'} + if strcmp(upper(gdat_data.gdat_params.equil),'FPG') + gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev; + gdat_params_eff.data_request = [{'FPG'},{'fax-bnd'},{'D3DD'}]; + gdat_data = gdat_d3d(gdat_data.shot,gdat_params_eff); + gdat_data.gdat_params = gdat_params_prev; + gdat_data.label = 'psi\_axis-psi\_edge'; + gdat_data.data_fullpath = [gdat_data.data_fullpath ' yields only psi\_axis-psi\_edge from FPG/fax-bnd']; + else + [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data); + % since Lpf depends on time, need to load all first and then loop over time for easier mapping + [psi_tree,e]=rdaD3D_eff(shot,DIAG,'PFL',exp_name_eff); + ndimrho = size(psi_tree.data,2); + if ndimrho==NTIME_Lpf + % data seems to be transposed + ndimrho = size(psi_tree.data,1); + itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t + else + itotransposeback = 0; + end + psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback); + ijnan=find(isnan(psi_tree.value)); + psi_tree.value(ijnan)=0; + gdat_data.dim{1} = gdat_data.t; + gdat_data.dimunits{1} = 'time [s]'; + for it=1:NTIME + Lpf1 = Lpf1_t(it); + psi_it=psi_tree.value(it,Lpf1:-1:1)'; + if strcmp(data_request_eff,'psi_axis') + gdat_data.data(it)= psi_it(1); + elseif strcmp(data_request_eff,'psi_edge') + gdat_data.data(it)= psi_it(end); + else + end + end + gdat_data.units = psi_tree.units; + gdat_data.data_fullpath = [DIAG '/PFL extract ' data_request_eff]; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'rhotor', 'rhotor_edge', 'rhotor_norm', 'rhovol', 'volume_rho'} + if strcmp(upper(gdat_data.gdat_params.equil),'FPG') + gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev; + gdat_params_eff.data_request = [{'FPG'},{'Vol'},{'D3DD'}]; + gdat_data = gdat_d3d(gdat_data.shot,gdat_params_eff); + gdat_data.gdat_params = gdat_params_prev; + gdat_data.label = 'Vol'; + gdat_data.data_fullpath = [gdat_data.data_fullpath ' yields only edge Volume from FPG/Vol']; + return + end + [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data); + % since Lpf depends on time, need to load all first and then loop over time for easier mapping + [phi_tree,e]=rdaD3D_eff(shot,DIAG,'TFLx',exp_name_eff); + ndimrho = size(phi_tree.data,2); + if ndimrho==NTIME_Lpf + % data seems to be transposed + ndimrho = size(phi_tree.data,1); + itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t + else + itotransposeback = 0; + end + phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback); + ijnan=find(isnan(phi_tree.value)); + phi_tree.value(ijnan)=0; + [Vol,e]=rdaD3D_eff(shot,DIAG,'Vol',exp_name_eff); + Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback); + [psi_tree,e]=rdaD3D_eff(shot,DIAG,'PFL',exp_name_eff); + psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback); + % + switch data_request_eff + case {'volume_rho', 'rhovol'} + for it=1:NTIME + Lpf1 = Lpf1_t(it); + psi_it=psi_tree.value(it,Lpf1:-1:1)'; + psi_axis(it)= psi_it(1); + psi_lcfs(it)= psi_it(end); + gdat_data.x(:,it) = sqrt(abs((psi_it-psi_axis(it)) ./(psi_lcfs(it)-psi_axis(it)))); + gdat_data.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)'; + % gdat_data.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi + gdat_data.rhovolnorm(:,it) = sqrt(abs(gdat_data.vol(:,it) ./ gdat_data.vol(end,it))); + end + gdat_data.dim{1} = gdat_data.x; + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits{1} = 'rhopolnorm'; + gdat_data.dimunits{2} = 'time [s]'; + if strcmp(data_request_eff,'volume_rho') + gdat_data.data = gdat_data.vol; + gdat_data.units = Vol.units; + else strcmp(data_request_eff,'rhovol') + gdat_data.data = gdat_data.rhovolnorm; + gdat_data.units = ' '; + end + case {'rhotor', 'rhotor_edge', 'rhotor_norm'} + b0=gdat(shot,'b0'); + gdat_data.b0 = interpos(b0.t,b0.data,gdat_data.t,-1); + for it=1:NTIME + Lpf1 = Lpf1_t(it); + gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)'; + gdat_data.rhotornorm(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.phi(end,it))); + gdat_data.rhotor(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.b0(it) ./ pi)); + psi_it=psi_tree.value(it,Lpf1:-1:1)'; + psi_axis(it)= psi_it(1); + psi_lcfs(it)= psi_it(end); + gdat_data.x(:,it) = sqrt(abs((psi_it-psi_axis(it)) ./(psi_lcfs(it)-psi_axis(it)))); + end + if strcmp(data_request_eff,'rhotor') + gdat_data.data = gdat_data.rhotor; + gdat_data.units = 'm'; + gdat_data.dim{1} = gdat_data.x; + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits{1} = 'rhopolnorm'; + gdat_data.dimunits{2} = 'time [s]'; + elseif strcmp(data_request_eff,'rhotor_edge') + gdat_data.data = gdat_data.rhotor(end,:); + gdat_data.units = 'm'; + gdat_data.dim{1} = gdat_data.t; + gdat_data.dimunits{1} = 'time [s]'; + elseif strcmp(data_request_eff,'rhotor_norm') + gdat_data.data = gdat_data.rhotornorm; + gdat_data.units = ' '; + gdat_data.dim{1} = gdat_data.x; + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits{1} = 'rhopolnorm'; + gdat_data.dimunits{2} = 'time [s]'; + end + end + gdat_data.data_fullpath = [DIAG '/PFL extract in .data: ' data_request_eff]; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'sxr'} + % sxr from B by default or else if 'source','else' is provided + if ~isfield(gdat_data.gdat_params,'freq')|| isempty(gdat_data.gdat_params.freq) + gdat_data.gdat_params.freq = 'slow'; + end + if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) + gdat_data.gdat_params.source = 'G'; + elseif length(gdat_data.gdat_params.source)>1 || ... + upper(gdat_data.gdat_params.source)<'F' || upper(gdat_data.gdat_params.source)>'M' + if gdat_data.gdat_params.nverbose>=1 + warning(['source = ' gdat_data.gdat_params.source ' not expected with data_request= ' data_request_eff]) + end + return + end + dim1_len = 'from_chord_nb'; % thus up to max(chords_ok_nb) + if ~isfield(gdat_data.gdat_params,'camera') || isempty(gdat_data.gdat_params.camera) + gdat_data.gdat_params.camera = []; + elseif ischar(gdat_data.gdat_params.camera) && strcmp(gdat_data.gdat_params.camera,'central') + gdat_data.gdat_params.source = 'J'; + gdat_data.gdat_params.camera = 49; + elseif isnumeric(gdat_data.gdat_params.camera) + % ok keep the array, but first dim to contain just the related chords + dim1_len='nb_of_chords'; + else + if gdat_data.gdat_params.nverbose>=1 + warning(['camera = ' gdat_data.gdat_params.camera ' not expected with data_request= ' data_request_eff]) + end + return + end + if length(gdat_data.gdat_params.camera)==1; dim1_len='nb_of_chords'; end + gdat_data.gdat_params.source = upper(gdat_data.gdat_params.source); + % + if ~isfield(gdat_data.gdat_params,'time_interval') + gdat_data.gdat_params.time_interval = []; + end + [aa,bb]=unix(['ssh ' '-o "StrictHostKeyChecking no" sxd3d20.d3d.ipp.mpg.de WhichSX ' num2str(shot) ' ' ... + upper(gdat_data.gdat_params.source)]); + chords_ok_diags=regexpi(bb,'(?<chord>[A-Z]_[0-9]+)\saddress:\s[0-9]+\sDiag:\s(?<diag>[A-Z]+)','names'); + if isempty(chords_ok_diags) + chords_ok_nb = []; + else + for i=1:length(chords_ok_diags) + chords_ok_nb(i) = str2num(chords_ok_diags(i).chord(3:end)); + end + end + if isempty(gdat_data.gdat_params.camera); + gdat_data.gdat_params.camera = chords_ok_nb; + else + for i=1:length(gdat_data.gdat_params.camera) + ij=find(chords_ok_nb==gdat_data.gdat_params.camera(i)); + if ~isempty(ij) + chords_ok_diags_tmp(i) = chords_ok_diags(ij); + else + % chord not in use + chords_ok_diags_tmp(i).chord = []; + chords_ok_diags_tmp(i).diag = []; + end + end + chords_ok_diags = chords_ok_diags_tmp; + chords_ok_nb = gdat_data.gdat_params.camera; + end + exp_name_eff = 'D3DD'; + icount = 0; + nnth = 1; + if isnumeric(gdat_data.gdat_params.freq) && gdat_data.gdat_params.freq>1; + nnth = floor(gdat_data.gdat_params.freq+0.5); + gdat_data.gdat_params.freq = nnth; + end + for i=1:length(gdat_data.gdat_params.camera) + if ischar(gdat_data.gdat_params.freq) && strcmp(gdat_data.gdat_params.freq,'slow'); chords_ok_diags(i).diag = 'SSX'; end + if ~isempty(chords_ok_diags(i).diag) && ~isempty(chords_ok_diags(i).chord) + [a,e]=rdaD3D_eff(shot,chords_ok_diags(i).diag,chords_ok_diags(i).chord,exp_name_eff,gdat_data.gdat_params.time_interval); + else + a.data = []; + a.t = []; + end + if ~isempty(a.data) + icount = icount + 1; + if icount == 1 + % first time has data + % assume all chords have same time base even if from different shotfile + % time missing one point + if length(a.value) == length(a.t)+1 + a.t=linspace(a.range(1),a.range(2),length(a.value)); + a_time.index(2) = length(a.value); + end + gdat_data.t = a.t(1:nnth:end); + gdat_data.units = a.units; + if strcmp(dim1_len,'from_chord_nb') + gdat_data.data = NaN*ones(max(chords_ok_nb),length(gdat_data.t)); + gdat_data.dim{1} = [1:max(chords_ok_nb)]; % simpler to have index corresponding to chord number, except for central + else + gdat_data.data = NaN*ones(length(chords_ok_nb),length(gdat_data.t)); + gdat_data.dim{1} = chords_ok_nb; + end + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits = [{'chord nb'}; {'s'}]; + gdat_data.data_fullpath = ['sxr from source=''' gdat_data.gdat_params.source ' uses shotfiles SX' ... + setdiff(unique(strvcat(chords_ok_diags(:).diag)),'SX')']; + gdat_data.label = ['SXR/' upper(gdat_data.gdat_params.source)]; + end + try + if strcmp(dim1_len,'from_chord_nb') + gdat_data.data(chords_ok_nb(i),:) = a.data(1:nnth:end); + else + gdat_data.data(i,:) = a.data(1:nnth:end); + end + catch + if (gdat_params.nverbose>=1); disp(['problem with ' chords_ok_diags(i).diag ' ; ' chords_ok_diags(i).chord]); end + + end + else + % add fields not yet defined in case all cases have empty data + end + end + gdat_data.chords = chords_ok_nb; + if length(gdat_data.dim)>=1 + gdat_data.x = gdat_data.dim{1}; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'transp'} + % most of the times the exp for the shotfile should be provided + shotfile_exp_eff = gdat_params.exp_name; + diagname='TRA'; + TRANSP_signals; + for i=1:size(transp_sig,1) + if strcmp(lower(transp_sig{i,2}),'signal') || strcmp(lower(transp_sig{i,2}),'signal-group') + try + eval(['[gdat_data.' transp_sig{i,1} ',e]=rdaD3D_eff(shot,diagname,''' transp_sig{i,1} ''',shotfile_exp_eff);']); + catch + eval(['gdat_data.' transp_sig{i,1} '=[];']); + end + elseif strcmp(lower(transp_sig{i,2}),'area-base') + clear adata_area + try + [adata_area]=sf2ab(diagname,shot,transp_sig{i,1},'-exp',shotfile_exp_eff); + catch + adata_area.value = cell(0); + end + eval(['gdat_data.' transp_sig{i,1} '=adata_area;']); + elseif strcmp(lower(transp_sig{i,2}),'time-base') + clear adata_time + try + [adata_time]=sf2tb(diagname,shot,transp_sig{i,1},'-exp',shotfile_exp_eff); + catch + adata_time.value = cell(0); + end + eval(['gdat_data.' transp_sig{i,1} '=adata_time;']); + end + end + % copy TIME to .t + if isfield(gdat_data,'TIME') && isfield(gdat_data.TIME,'value') + gdat_data.t = gdat_data.TIME.value; + gdat_data.dim{1} = gdat_data.t, + gdat_data.dimunits{1} = gdat_data.TIME.unit; + end + + otherwise + if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_d3d']); end + error_status=901; + return + end + +else + if (gdat_params.nverbose>=1); warning(['D3D method=' mapping_for_D3D.method ' not known yet, contact Olivier.Sauter@epfl.ch']); end + error_status=602; + return +end + +gdat_data.gdat_params.help = d3d_help_parameters(fieldnames(gdat_data.gdat_params)); +gdat_data.mapping_for.d3d = mapping_for_d3d; +gdat_params = gdat_data.gdat_params; +error_status=0; + +return + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% functions for portions called several times + +function [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data); +% get basic params to be able to read results in EQ-like shotfiles +% M_Rmesh,N_Zmesh only needed for equil when 2D quantities are required +% +shot=gdat_data.shot; +exp_name_eff = gdat_data.gdat_params.exp_name; +gdat_data.gdat_params.equil = upper(gdat_data.gdat_params.equil); +DIAG = gdat_data.gdat_params.equil; % 'EQI' by default at this stage, should become EQH? +M_Rmesh_par = sf2par(DIAG,shot,'M','PARMV'); +M_Rmesh = M_Rmesh_par.value + 1; % nb of points +N_Zmesh_par = sf2par(DIAG,shot,'N','PARMV'); +N_Zmesh = N_Zmesh_par.value + 1; % nb of points +NTIME_par = sf2par(DIAG,shot,'NTIME','PARMV'); +NTIME = NTIME_par.value; % nb of points +Lpf_par = rdaD3D_eff(shot,DIAG,'Lpf',exp_name_eff); +% since June, nb of time points in EQ results is not consistent with NTIME and time +% It seems the first NTIME points are correct, so use this explicitely +NTIME_Lpf = length(Lpf_par.value); +if (NTIME < NTIME_Lpf) + if gdat_data.gdat_params.nverbose>=3; disp('WARNING: nb of times points smaller then equil results, use first NTIME points'); end +elseif (NTIME > NTIME_Lpf) + if gdat_data.gdat_params.nverbose >= 1 + disp('ERROR: nb of times points LARGER then equil results') + disp('this is unexpected, so stop there and ask Olivier.Sauter@epfl.ch') + end + return +end +Lpf_tot = Lpf_par.value(1:NTIME); % nb of points: 100000*nb_SOL points + nb_core +Lpf_SOL = fix(Lpf_tot/100000); +Lpf1_t = mod(Lpf_tot,100000)+1; % nb of Lpf points + +[equil_time,e]=rdaD3D_eff(shot,DIAG,'time',exp_name_eff); +gdat_data.t = equil_time.value(1:NTIME); diff --git a/crpptbx/D3D/get_signal_d3d.m b/crpptbx/D3D/get_signal_d3d.m new file mode 100644 index 00000000..e1e1a7eb --- /dev/null +++ b/crpptbx/D3D/get_signal_d3d.m @@ -0,0 +1,28 @@ +function [gdat_data,error_status] = get_signal_d3d(signal,gdat_data); + +gdat_data.data = mdsvalue(signal); +error_status = 1; + +nbdims=numel(setdiff(size(gdat_data.data),1)); +if ~isempty(gdat_data.data) && ~ischar(gdat_data.data) + gdat_data.units = mdsvalue(['units_of(' signal ')']); + for i=1:nbdims + gdat_data.dim{i} = mdsvalue(['dim_of(' signal ',' num2str(i-1) ')']); + % gdat_data.dimunits{i} = mdsvalue(['dimunits_of(' signal ',' num2str(i-1) ')']); + gdat_data.dimunits{i} = 'to get'; + end + if nbdims==1 + gdat_data.t = gdat_data.dim{1}; + gdat_data.x = []; + else + gdat_data.t = gdat_data.dim{1}; + gdat_data.x = gdat_data.dim{2}; + end + error_status = 0; +else + gdat_data.dim = []; + gdat_data.dimunits = []; + gdat_data.t = []; + gdat_data.x = []; +end + diff --git a/crpptbx/gdat.m b/crpptbx/gdat.m index b80bb3eb..a02c2f62 100644 --- a/crpptbx/gdat.m +++ b/crpptbx/gdat.m @@ -87,6 +87,8 @@ else default_machine = 'tcv'; elseif ~isempty(regexpi(hostname,'rzg.mpg')) default_machine = 'aug'; + elseif ~isempty(regexpi(hostname,'gat.com')) + default_machine = 'd3d'; end end machine_eff = default_machine; diff --git a/crpptbx/gdat_plot.m b/crpptbx/gdat_plot.m index 60183c72..5f12e909 100644 --- a/crpptbx/gdat_plot.m +++ b/crpptbx/gdat_plot.m @@ -76,7 +76,7 @@ if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty( hylabel=ylabel(ylabel_eff); end if iscell(gdat_data.label) && length(gdat_data.label)>=2; - legend(gdat_data.label) + legend(gdat_data.label,'interpreter','none'); end zoom on; end diff --git a/crpptbx/gdatpaths.m b/crpptbx/gdatpaths.m index ca08087f..86be4fac 100644 --- a/crpptbx/gdatpaths.m +++ b/crpptbx/gdatpaths.m @@ -9,6 +9,7 @@ a=a(1:ii(end)); machines=[{'JET'} {'TCV'} {'AUG'} {'D3D'}]; machines=[{'JET'} {'TCV'} {'AUG'} {'KSTAR'}]; +machines=[{'JET'} {'TCV'} {'AUG'} {'D3D'} {'KSTAR'}]; for i=1:length(machines) addpath([a machines{i}]) diff --git a/crpptbx/get_data_request_names_from_gdat_xxx.m b/crpptbx/get_data_request_names_from_gdat_xxx.m index ad3fe967..3c35af49 100644 --- a/crpptbx/get_data_request_names_from_gdat_xxx.m +++ b/crpptbx/get_data_request_names_from_gdat_xxx.m @@ -7,7 +7,7 @@ function [data_request_names] = get_data_request_names_from_gdat_xxx(machine_in) data_request_names = []; -expected_machines = [{'aug'}, {'jet'}, {'tcv'}]; % substrutures created for these at this stage (all means all of these) +expected_machines = [{'aug'}, {'d3d'}, {'jet'}, {'tcv'}]; % substrutures created for these at this stage (all means all of these) for j=1:length(expected_machines) data_request_names.(expected_machines{j}) = []; end -- GitLab