diff --git a/crpptbx/TCV/gdat_tcv.m b/crpptbx/TCV/gdat_tcv.m index b45ace512c71b7385aca2d48c0ef1beb013ef9c2..6ce1bf373538ae622e93f793a3a8c2d93fb1ba00 100644 --- a/crpptbx/TCV/gdat_tcv.m +++ b/crpptbx/TCV/gdat_tcv.m @@ -2,7 +2,7 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(shot,data_req % % 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. +% 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 @@ -48,7 +48,7 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(shot,data_req % % Examples: % (should add working examples for various machines (provides working shot numbers for each machine...)) -% +% % [a1,a2]=gdat; % a2.data_request = 'Ip'; % a3=gdat(48836,a2); % gives input parameters as a structure, allows to call the same for many shots @@ -109,7 +109,7 @@ gdat_data.data_fullpath = []; % Treat inputs: ivarargin_first_char = 3; data_request_eff = ''; -if nargin>=2 && ischar(data_request); +if nargin>=2 && ischar(data_request); data_request_eff = data_request; data_request = data_request; % should not lower until see if keyword end @@ -129,7 +129,8 @@ do_mdsopen_mdsclose = 1; if nargin>=1 if isempty(shot) % means mdsopen(shot) already performed - shot = mdsipmex(2,'$SHOT'); + shot_mds = mdsipmex(2,'$SHOT'); + if isnumeric(shot_mds); shot = shot_mds; end gdat_data.shot = shot; do_mdsopen_mdsclose = 0; if ischar(shot) || isempty(shot) @@ -142,7 +143,7 @@ if nargin>=1 warning(['shot cannot be opened']); end end - return + if ~strcmp(data_request,'ids'); return; end % empty shot return empty ids so valid command end elseif isnumeric(shot) gdat_data.shot = shot; @@ -234,7 +235,7 @@ if ischar(data_request_eff) || length(data_request_eff)==1 else ij=[]; end -if ~isempty(ij); +if ~isempty(ij); gdat_data.gdat_request = data_request_names_all{ij}; gdat_params.data_request = gdat_data.gdat_request; if isfield(data_request_names.all.(data_request_names_all{ij}),'description') && ~isempty(data_request_names.all.(data_request_names_all{ij}).description) @@ -281,10 +282,10 @@ end if liuqe_version_eff==2 || liuqe_version_eff==3 substr_liuqe_tcv_eq = num2str(liuqe_version_eff); end - + % special treatment for model shot=-1 or preparation shot >=100'000 begstr = ''; -if shot==-1 || (shot>=100000 && shot < 200000) || liuqe_version==-1 +if ~isempty(shot) && (shot==-1 || (shot>=100000 && shot < 200000) || liuqe_version==-1 ) % requires FBTE liuqe_version_eff = -1; begstr = 'tcv_eq( "'; @@ -453,7 +454,7 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi') gdat_data.data = aatmp.data; gdat_data.dim = aatmp.dim; nbdims = length(gdat_data.dim); - if mapping_for_tcv.timedim==-1; + if mapping_for_tcv.timedim==-1; % try to find time dim from units idim_non1 = []; len_dim = []; for i=1:length(aatmp.dimunits) @@ -492,7 +493,7 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi') gdat_data.dimunits = aatmp.dimunits; if mapping_for_tcv.gdat_timedim>0 && mapping_for_tcv.gdat_timedim ~= mapping_for_tcv.timedim % shift timedim to gdat_timedim data(i,j,...itime,k,...) -> data(i,inewtime,j,...,k,...) - % note that this means that gdat_data.x and gdat_data.t are same and correct, + % note that this means that gdat_data.x and gdat_data.t are same and correct, % only .data, .dim and .dimunits need to be changed iprev=[1:nbdims]; ij=find(dim_nontim>mapping_for_tcv.gdat_timedim-1); @@ -544,7 +545,7 @@ elseif strcmp(mapping_for_tcv.method,'expression') ijdim=find(strcmp(tmp_fieldnames,'dim')==1); if ~isempty(ijdim) nbdims = length(gdat_data.dim); - if mapping_for_tcv.timedim==-1; + if mapping_for_tcv.timedim==-1; mapping_for_tcv.timedim = nbdims; if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_tcv.timedim = nbdims-1; end end @@ -581,8 +582,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'a_minor','rgeom','r_geom','a_minor_rho','r_geom_rho','rgeom_rho'} - if strcmp(data_request_eff,'r_geom'); data_request_eff = 'rgeom'; end - if strcmp(data_request_eff,'r_geom_rho'); data_request_eff = 'rgeom_rho'; end + if strcmp(data_request_eff,'r_geom'); data_request_eff = 'rgeom'; end + if strcmp(data_request_eff,'r_geom_rho'); data_request_eff = 'rgeom_rho'; end % compute average minor or major radius (on z=zaxis normally) if liuqe_matlab==0 nodenameeff=['tcv_eq(''r_max_psi'',''LIUQE' substr_liuqe_tcv_eq ''')']; @@ -593,7 +594,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end return - end + end nodenameeff2=['tcv_eq(''r_min_psi'',''LIUQE' substr_liuqe_tcv_eq ''')']; rminpsi=tdi(nodenameeff2); ijnan = find(isnan(rminpsi.data)); @@ -669,15 +670,15 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.data_fullpath=[nodenameeff ' and ' nodenameeff2]; else gdat_data.data = aatmp.data(end,:)'; - gdat_data.dim = aatmp.dim(2); - aatmp.dimunits = aatmp.dimunits(2); + gdat_data.dim = aatmp.dim(2); + aatmp.dimunits = aatmp.dimunits(2); gdat_data.data_fullpath=[nodenameeff]; end gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim}; gdat_data.units = aatmp.units(end); gdat_data.dimunits = aatmp.dimunits; end - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'zgeom','z_geom'} % compute average minor or major radius (on z=zaxis normally) @@ -697,7 +698,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end return - end + end if strcmp(data_request_eff,'zgeom') gdat_data.data=0.5.*(max(zcontour.data,[],1) + min(zcontour.data,[],1)); gdat_data.data_fullpath=['(max+min)/2 of ' nodenameeff]; @@ -711,7 +712,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if any(strcmp(fieldnames(zcontour),'units')) gdat_data.units = zcontour.units; end - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'b0'} % B0 at R0=0.88 @@ -721,7 +722,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') % expect liuqe or liuqe.m to have liuqe time-base, otherwise give full time base else gdat_data.gdat_params.source = 'iphi'; - end + end if liuqe_version_eff==-1 nodenameeff = 'tcv_eq("BZERO","FBTE")'; tracetdi=tdi(nodenameeff); @@ -747,7 +748,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') return end gdat_data.data_fullpath=[nodenameeff]; - gdat_data.dim = tracetdi.dim; + gdat_data.dim = tracetdi.dim; gdat_data.t = gdat_data.dim{1}; if any(strcmp(fieldnames(tracetdi),'units')) gdat_data.units = tracetdi.units; @@ -755,14 +756,14 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.dimunits = tracetdi.dimunits; gdat_data.request_description = ['vacuum magnetic field at R0=' num2str(r0exp) 'm; COCOS=17']; gdat_data.r0 = r0exp; - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'betan'} % 100*beta / |Ip[MA] * B0[T]| * a[m] % get B0 from gdat_tcv, without re-opening the shot and using the same parameters except data_request % easily done thanks to structure call for options params_eff = gdat_data.gdat_params; - params_eff.data_request='b0'; + params_eff.data_request='b0'; b0=gdat_tcv([],params_eff); % note: no need to set .doplot=0 since gdat_tcv does not call gdat_plot in any case params_eff.data_request='ip'; ip=gdat_tcv([],params_eff); @@ -788,16 +789,16 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.data_fullpath='100*beta/ip*1e6*b0*a_minor, each from gdat_tcv'; gdat_data.units = ''; gdat_data.dimunits{1} = 's'; - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'cxrs','cxrs_rho'} %not yet finished, just started % load typical data from cxrs, Ti, ni, vtori and vpoli (if available), as well as zeff from cxrs % if 'fit' option is added: 'fit',1, then the fitted profiles are returned - % + % % sub_nodes names from CXRS_get_profiles function, lower case when used in gdat sub_nodes = {'Ti','vTor','vPol','nC','Zeff'}; % first node is also copied into data, choose "default' one - sub_nodes_out = lower({'Ti','vTor','vPol','ni','Zeff'}); % + sub_nodes_out = lower({'Ti','vTor','vPol','ni','Zeff'}); % sub_nodes_units = {'eV','m/s','m/s','m^{-3}',''}; % first node is also copied into data, choose "default' one % use A. Karpushov routine to get profiles and then copy the data or the fitted profiles aa=CXRS_get_profiles; cxrs_params = aa.param; @@ -847,7 +848,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') cxrs_params.prof.Ti.taus = fit_tension.ti; cxrs_params.prof.vi.taus = fit_tension.vi; cxrs_params.prof.nc.taus = fit_tension.ni; - cxrs_params.prof.zeff.taus = fit_tension.zeff; + cxrs_params.prof.zeff.taus = fit_tension.zeff; cxrs_params.k_plot = cxrs_plot; cxrs_profiles = CXRS_get_profiles(shot,source,time_interval,cxrs_params); inb_times = length(cxrs_profiles.Times); @@ -905,7 +906,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.(sub_eff_out).raw.error_bar_rho(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.dx; gdat_data.(sub_eff_out).raw.cxrs_system(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.sys; gdat_data.time_interval{it} = cxrs_profiles.(sub_eff){it}.t_lim; - end + end end gdat_data.(sub_eff_out).units = sub_nodes_units{i}; if i==1 @@ -1023,7 +1024,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.gdat_params.cocos = cocos_out; end [eqdsk_cocosout, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdsk_cocos_in,[cocos_in cocos_out]); - % for several times, use array of structure for eqdsks, + % for several times, use array of structure for eqdsks, % cannot use it for psi(R,Z) in .data and .x since R, Z might be different at different times, % so project psi(R,Z) on Rmesh, Zmesh of 1st time if length(time) > 1 @@ -1069,7 +1070,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') 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 {'halpha'} channels = -1; @@ -1084,28 +1085,51 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if isfield(params_eff,'source') && ~isempty(params_eff.source) ids_top_name = params_eff.source; else - disp('need an ids name in ''source'' parameter') - return - end - try - equil_empty=ids_gen(ids_top_name); - catch - warning(['cannot create empty ids: ids_gen(''' ids_top_name ''')']); - return + ids_top_name = []; + disp('need an ids name in ''source'' parameter'); + disp('check substructure gdat_params.source_available for an ids list'); + end + ids_gen_ok = exist('ids_gen'); + equil_empty = struct([]); + if ids_gen_ok ~= 2 + ids_struct_saved = 'ids_structures_20190312.mat'; + if ~exist(ids_struct_saved,'file') + warning(['function ids_gen not available neither file ids_structures_20190312.mat thus cannot create empty ids: ids_gen(''' ids_top_name ''')']); + return + else + load ids_structures_20190312.mat + if isfield(ids_structures,ids_top_name) + equil_empty = ids_structures.(ids_top_name); + else + if ~isempty(ids_top_name); + warning(['ids ''' ids_top_name ''' does not exist in ids_structures saved in ' ids_struct_saved]); + end + gdat_data.gdat_params.source_available = ids_list; + return + end + end + else + equil_empty = ids_gen(ids_top_name); end try - eval(['[ids_top,ids_top_description]=tcv_get_ids_' ids_top_name '(shot,equil_empty);']) - gdat_data.(ids_top_name) = ids_top; - gdat_data.([ids_top_name '_description']) = ids_top_description; - catch + if ~isempty(shot) + eval(['[ids_top,ids_top_description]=tcv_get_ids_' ids_top_name '(shot,equil_empty);']) + gdat_data.(ids_top_name) = ids_top; + gdat_data.([ids_top_name '_description']) = ids_top_description; + else + gdat_data.(ids_top_name) = equil_empty; + gdat_data.([ids_top_name '_description']) = ['shot empty so return default empty structure for ' ids_top_name]; + end + catch ME_tcv_get_ids + getReport(ME_tcv_get_ids) disp(['there is a problem with: tcv_get_ids_' ids_top_name ... ' , may be check if it exists in your path or test it by itself']) - gdat_data.(ids_top_name) = struct([]); - gdat_data.([ids_top_name '_description']) = struct([]); + gdat_data.(ids_top_name) = equil_empty; + gdat_data.([ids_top_name '_description']) = getReport(ME_tcv_get_ids); end - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - case {'icds'} + case {'pdens','icds'} % note: same time array for all at main.data level, then individual at .ec, .nbi levels % At this stage fill just eccd, later add nbi sources_avail = {'ec','nbi'}; % can be set in parameter source @@ -1172,14 +1196,15 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') % EC if strcmp(lower(source_icd.ec),'toray') if isempty(gdat_data.gdat_params.trialindx) - [ptot,icdtot] = astra_tcv_EC_exp(shot); % centralized function for toray nodes + % [ p_gyro,icdtot,ECH,ECCD,rho_dep,drho,pmax,icdmax,ECCD_w2] = astra_tcv_EC_exp(shot,time,cd_sign,smooth_pgyro,slices,toray_trialindex) + [pabs_gyro,icdtot,ECH,ECCD,rho_dep_pow,drho_pow,pmax,icdmax,ECCD_w2,rho_dep_icd,drho_icd] = astra_tcv_EC_exp(shot); % centralized function for toray nodes else - [ptot,icdtot] = astra_tcv_EC_exp(shot,[],[],[],[],[],gdat_data.gdat_params.trialindx); % centralized function for toray nodes + [pabs_gyro,icdtot,ECH,ECCD,rho_dep_pow,drho_pow,pmax,icdmax,ECCD_w2,rho_dep_icd,drho_icd] = astra_tcv_EC_exp(shot,[],[],[],[],[],gdat_data.gdat_params.trialindx); % centralized function for toray nodes end - data_fullpath = 'from toray nodes using astra_tcv_EC_exp(shot)'; + data_fullpath = ['from toray nodes using astra_tcv_EC_exp(shot), total icd per gyrotron including effective absorbed power, jcd etc in .extra']; ec_help = 'from toray icdint with extracting of effective Icd for given launcher depending on nb rays used'; if isfield(icdtot,'tgrid'); icdtot.t = icdtot.tgrid; end - if isfield(icdtot,'data'); + if isfield(icdtot,'data'); icdtot.data = icdtot.data * 1e6; % MA -> A icdtot.data = icdtot.data'; icdtot.data(:,end+1) = sum(icdtot.data,2); @@ -1249,7 +1274,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') % used in gdat_plot for spectrogram plot else gdat_data.gdat_params.nfft = 1024; - end + end % load n=1, 2 and 3 Bdot from magnetic measurements n3.data = []; if shot< 50926 @@ -1262,7 +1287,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') % gdat_data.gdat_params.source; else gdat_data.gdat_params.source = '23'; - end + end if length(gdat_data.gdat_params.source)>=2 && strcmp(gdat_data.gdat_params.source(1:2),'23') aaLFSz23_sect3=tdi('\atlas::DT196_MHD_001:channel_067'); aaLFSz23_sect11=tdi('\atlas::DT196_MHD_001:channel_075'); @@ -1302,7 +1327,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.n2_HFS.data = gdat_data.n2_HFS.data + aaHFSz0_sect11.data; gdat_data.data_fullpath=['\atlas::DT196_MHD_001:channel_083 -+ \atlas::DT196_MHD_001:channel_091 for n=1,2, LFS_sect3/11, z=0cm; _HFS' ... ' same for HFS: MHD_002:channel_034-+MHD_002:channel_038']; - end + end else disp('should not be here in ''mhd'', ask O. Sauter') return @@ -1314,7 +1339,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if length(n3.data)==length(n1.data); gdat_data.data(:,3) = reshape(n3.data,length(n3.data),1); end gdat_data.dim{1} = n1.dim{1}; gdat_data.t = gdat_data.dim{1}; - gdat_data.dim{2} = [1; 2; 3]; + gdat_data.dim{2} = [1; 2; 3]; gdat_data.dimunits{1} = n1.dimunits{1}; gdat_data.dimunits{2} = 'n number'; if shot>= 50926 @@ -1328,7 +1353,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.label = {'n odd','n even'}; % can be used in legend(gdat_data.label) end end - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'ne','te'} % ne or Te from Thomson data on raw z mesh vs (z,t) @@ -1337,7 +1362,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.gdat_params.edge = 0; end [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.edge,gdat_params.nverbose); - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'ne_rho', 'te_rho', 'nete_rho'} % if nete_rho, do first ne, then Te later (so fir stuff already done) @@ -1430,7 +1455,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.x=rho; % add grids_1d to have rhotor, etc if cxrs_rho was asked for gdat_data = get_grids_1d(gdat_data,2,1,gdat_params.nverbose); - + %%%%%%%%%%% % if nete_rho, copy data as .ne, get .te and put pe=e ne Te in data: if strcmp(data_request_eff(1:4),'nete') @@ -1799,7 +1824,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end return - end + end gdat_data.data = tracetdi.data; gdat_data.dim = tracetdi.dim; gdat_data.t = gdat_data.dim{2}; @@ -1808,7 +1833,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') rhopol_eff = ones(size(tracetdi.dim{1})); rhopol_eff(:) = sqrt(linspace(0,1,length(tracetdi.dim{1}))); gdat_data.dim{1} = rhopol_eff; - end + end gdat_data.x = gdat_data.dim{1}; gdat_data.dimunits{1} = 'rho_pol~sqrt(\psi_norm)'; gdat_data.dimunits{2} = 's'; @@ -1835,7 +1860,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end return - end + end gdat_data.data = tracetdi.data; gdat_data.dim = tracetdi.dim; gdat_data.t = gdat_data.dim{2}; @@ -1844,7 +1869,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') rhopol_eff = ones(size(tracetdi.dim{1})); rhopol_eff(:) = sqrt(linspace(0,1,length(tracetdi.dim{1}))); gdat_data.dim{1} = rhopol_eff; - end + end gdat_data.x = gdat_data.dim{1}; gdat_data.dimunits{1} = 'rho_pol~sqrt(\psi_norm)'; gdat_data.dimunits{2} = 's'; @@ -1884,7 +1909,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') [qfit,~,~,phi]=interpos(-43,q_rho.dim{1}(ij).^2,q_rho.data(ij,it),phi_tor.dim{1}.^2); % Phi=sigma_bp*sigma_rhothetaphi*(2pi)^(1-eBp)*int(q dpsi), dpsi=dpsi_norm * (psiedge-psaxis) % cocos=17: sigma_Bp=-1,sigma_rhothetaphi=1, eBp=1, (psiedge-psaxis)=-psiaxis - phi = -1.* phi .* (-psi_axis.data(it)); + phi = -1.* phi .* (-psi_axis.data(it)); phi_tor.data(:,it) = phi'; end end @@ -1913,7 +1938,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ij=find(isfinite(q_rho.data(:,it))); if ~isempty(ij) && numel(ij)>5 [qfit,~,~,phi]=interpos(-43,q_rho.dim{1}(ij).^2,q_rho.data(ij,it),phi_tor.dim{1}.^2); - phi = -1.* phi .* (-psi_axis.data(it)); + phi = -1.* phi .* (-psi_axis.data(it)); if only_edge phi_tor.data(end,it) = phi(end); else @@ -1929,7 +1954,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.units = phi_tor.units; if (length(gdat_data.dim)>=mapping_for_tcv.gdat_timedim); gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim}; end if (length(gdat_data.dim)>0); gdat_data.x = gdat_data.dim{1}; end - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'pprime', 'pprime_rho'} if liuqe_matlab==0 @@ -1990,7 +2015,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end return - end + end gdat_data.data = tracetdi.data.*0; gdat_data.dim = tracetdi.dim; gdat_data.t = gdat_data.dim{1}; @@ -2110,10 +2135,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.rhotor_edge = gdat_data.rhotor_edge'; else params_eff = gdat_data.gdat_params; - params_eff.data_request='phi_tor'; + params_eff.data_request='phi_tor'; phi_tor=gdat_tcv([],params_eff); params_eff = gdat_data.gdat_params; - params_eff.data_request='b0'; + params_eff.data_request='b0'; b0=gdat_tcv([],params_eff); gdat_data.t = phi_tor.t; ij=find(isfinite(b0.data)); @@ -2150,7 +2175,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.rhotor_edge = []; end end - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'rhovol','rho_vol','volume_rho','volume'} if strcmp(data_request_eff,'rho_vol'); data_request_eff = 'rhovol'; end @@ -2218,11 +2243,11 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end end gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim}; - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'rtc'} % load all real-time memory signals for various nodes - + %Get the data from mds and fill the data structure defined by %define_simulink_signals @@ -2277,7 +2302,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') else gdat_data.gdat_params.max_adcs = []; end - + % Get data for the source requested for isource=1:numel(gdat_data.gdat_params.source) switch gdat_data.gdat_params.source{isource} @@ -2288,13 +2313,13 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') mdsclose; mdsdisconnect; mdsopen(shot); - + case 'defined' DS = take_all_signals(shot); %Read from mds SDS_DS = define_simulink_signals(DS); %Put them in predifined structure % Add the .data and .t structure % there was a conflict with r 7471 - for ii=1:numel(SDS_DS) %iter over node + for ii=1:numel(SDS_DS) %iter over node node = ii; for jj=1:numel(SDS_DS{ii}) %iter over threads thread = jj; @@ -2321,7 +2346,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') tmp = tdi(expression); if isnumeric(tmp.data) && ~isempty(tmp.data) SDS_DS{ii}{jj}.(fieldnameslist{kk}).data =[SDS_DS{ii}{jj}.(fieldnameslist{kk}).data; tmp.data']; - + SDS_DS{ii}{jj}.(fieldnameslist{kk}).t = tmp.dim{1}; else fprintf('Warning node: %d thread: %d signal: %s ind %d not available\n', ii,jj,fieldnameslist{kk}, zz ); @@ -2334,7 +2359,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.rtc_defined = SDS_DS; mdsclose; mdsdisconnect; - + case 'all' mdsconnect('scd'); mdsopen('rtc', shot); @@ -2348,7 +2373,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') 8,4,256]; % Default signals initialization AS = init_all_signals(global_node_thread_signals); - % Put signals in standard data strucure (SDS) + % Put signals in standard data strucure (SDS) SDS_AS = define_simulink_signals(AS); % Add the .data and .t structure for ii=1:numel(SDS_AS) %iter over node @@ -2394,7 +2419,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') %to be added in case end end - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'sxr', 'mpx'} @@ -2485,7 +2510,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end case 'xtomo' - % so far allow string and array as 'camera' choices: + % so far allow string and array as 'camera' choices: % camera = [] (default, thus get XTOMOGetData defaults), 'central', [3 6] (camera numbers) camera_xtomo = []; channel_xtomo = []; @@ -2510,7 +2535,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if gdat_data.gdat_params.nverbose>=3; disp(['camera = ' gdat_data.gdat_params.camera ' not implemented']); end gdat_data.gdat_params.camera = []; end - + try if isempty(gdat_data.gdat_params.time_interval); [sig,t,Sat_channel,cat_default] = XTOMOGetData(shot,0,4,camera_xtomo,[]); @@ -2548,9 +2573,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') warning(['camera = ' gdat_data.gdat_params.camera ' not expected with data_request= ' data_request_eff]) end return - + end - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'ttprime', 'ttprime_rho'} @@ -2579,7 +2604,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'profnerho','profterho'} % for backward compatibility but corresponds to ne_rho with param.fit_type='auto' (TCV special) - % + % nodenameeff=['\results::THOMSON.PROFILES.AUTO:' data_request_eff(5:6)]; nodenameeff_vers = [nodenameeff ':version_num']; avers = tdi(nodenameeff_vers); @@ -2633,7 +2658,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') error_status=901; return end - + else if (gdat_params.nverbose>=1); warning(['TCV method=' mapping_for_tcv.method ' not known yet, contact Olivier.Sauter@epfl.ch']); end error_status=602; @@ -2706,7 +2731,7 @@ elseif strcmp(maintracename_eff,'thomson_edge') else if (nverbose>=1); disp('bad input in get_fir_thom_rat_data'); end return -end +end if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) && any(isfinite(tracefirrat.data)) ... && ~isempty(tracefirrat.dim) && ~isempty(tracefirrat.dim{1}) @@ -2726,7 +2751,7 @@ catch gdat_data.data_raw = gdat_data.data; gdat_data.error_bar_raw = gdat_data.error_bar; end - + if (nverbose>=1) && shot<100000 warning('Problems with \results::thomson:times') warning(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff]) diff --git a/crpptbx/TCV/tcv_requests_mapping.m b/crpptbx/TCV/tcv_requests_mapping.m index bb49d28ad0655d1eb61a391e8e9476a56f530e0b..c4f6821c94a8d3923e92db0fa0a6194c1a2f8e28 100644 --- a/crpptbx/TCV/tcv_requests_mapping.m +++ b/crpptbx/TCV/tcv_requests_mapping.m @@ -229,6 +229,10 @@ switch lower(data_request) 'gdat_tmp2=gdat_tcv(shot,params_eff);ij=find(gdat_tmp2.data==0);gdat_tmp2.data(ij)=NaN;' ... 'tmp_data2=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ... 'gdat_tmp.data = gdat_tmp.data./(tmp_data2+1e-5);']; + case {'pdens'} + mapping.timedim = 1; + mapping.label = 'various Pdens'; + mapping.method = 'switchcase'; case {'phi_tor', 'phitor', 'toroidal_flux'} mapping.timedim = 2; mapping.method = 'tdiliuqe'; @@ -483,7 +487,7 @@ switch lower(data_request) 'gdat_tmp=gdat_tcv([],params_eff); ' ... 'params_eff = gdat_data.gdat_params;params_eff.data_request=''\magnetics::iplasma:trapeze''; ' ... 'aa=gdat_tcv([],params_eff);it=find(abs(aa.data)<10e3);it2=iround_os(gdat_tmp.t,aa.t(it));gdat_tmp.data(it2)=NaN;']; - + % $$$ case '' % $$$ mapping.timedim = 1; % $$$ mapping.label = data_request; @@ -493,7 +497,7 @@ switch lower(data_request) mapping.label = data_request; mapping.method = 'tdi'; % 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)