diff --git a/matlab/JET/jet_requests_mapping.m b/matlab/JET/jet_requests_mapping.m index 4ed7576f49814d24d16e7af2c3f8655fe71c2749..664bf4be1a85890b5eb2aebb0e9b888cb37ad249 100644 --- a/matlab/JET/jet_requests_mapping.m +++ b/matlab/JET/jet_requests_mapping.m @@ -340,10 +340,20 @@ switch lower(data_request) mapping.timedim = 1; mapping.method = 'signal'; mapping.expression = [{'PPF'},{'EFIT'},{'RMAG'}]; + case {'rnt', 'neutron_rate'} + mapping.label = 'neutron\_rate'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'PPF'},{'TIN'},{'RNT'}]; case {'sxr', 'sxr_s40', 'sxr_htv'} mapping.timedim = 1; mapping.gdat_timedim = 2; mapping.method = 'switchcase'; + case {'tbeo'} + mapping.label = 'TBEO'; + mapping.timedim = 1; + mapping.method = 'signal'; + mapping.expression = [{'PPF'},{'EDG8'},{'TBEO'}]; case 'te' mapping.timedim = 2; mapping.label = 'Te'; diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m index 11023819e8e91a2c040cf5ba609dc47a6a48855f..41907f55ceb05140887dd8ef60ba518361cf8774 100644 --- a/matlab/TCV/gdat_tcv.m +++ b/matlab/TCV/gdat_tcv.m @@ -312,10 +312,10 @@ end mapping_for_tcv = tcv_requests_mapping(data_request_eff,shot); gdat_data.label = mapping_for_tcv.label; -% special treatment for model shot=-1 or preparation shot >=100'000 +% special treatment for model shot=-1 or preparation shot >=900'000 begstr = ''; if (iscell(mapping_for_tcv.expression) || isempty(strfind(mapping_for_tcv.expression,'\rtc::'))) && ... - ~isempty(shot) && (shot==-1 || (shot>=100000 && shot < 200000) || liuqe_version==-1 ) + ~isempty(shot) && (shot==-1 || (shot>=900000 && shot <= 999999 ) || liuqe_version==-1 ) % requires FBTE liuqe_version_eff = -1; begstr = 'tcv_eq( "'; @@ -829,6 +829,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') case {'b0'} % B0 at R0=0.88 r0exp=0.88; + added_correction_str=''; if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source) ... && length(gdat_data.gdat_params.source)>=5 && strcmp(lower(gdat_data.gdat_params.source(1:5)),'liuqe') % expect liuqe or liuqe.m to have liuqe time-base, otherwise give full time base @@ -849,7 +850,6 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') else if liuqe_matlab==0 added_correction = 1.0; % correction already in liuqe.f - added_correction_str = ['']; nodenameeff = ['tcv_eq(''BZERO'',''LIUQE' substr_liuqe_tcv_eq ''')']; else if isempty(substr_liuqe); substr_liuqe = '_1'; end @@ -1453,6 +1453,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') rethrow(ME); end end + % set default inputs if not given if ~isfield(gdat_data.gdat_params,'error_bar') || isempty(gdat_data.gdat_params.error_bar) gdat_data.gdat_params.error_bar = 'delta'; end @@ -1586,9 +1587,17 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') % EC if strcmp(lower(source_icd.ec),'toray') if isempty(gdat_data.gdat_params.trialindx) - [pabs_gyro,icdtot,pow_dens,currentdrive_dens,rho_dep_pow,drho_pow,pmax,icdmax,currentdrive_dens_w2,rho_dep_icd,drho_icd] = astra_tcv_EC_exp(shot); % centralized function for toray nodes + % centralized function for toray nodes + [pabs_gyro,icdtot,pow_dens,currentdrive_dens,rho_dep_pow,drho_pow,... + pmax,icdmax,currentdrive_dens_w2,rho_dep_icd,drho_icd,... + ~,power_integrated,currentdrive_integrated] = ... + astra_tcv_EC_exp(shot,[],[],[],[],[],[],1); else - [pabs_gyro,icdtot,pow_dens,currentdrive_dens,rho_dep_pow,drho_pow,pmax,icdmax,currentdrive_dens_w2,rho_dep_icd,drho_icd] = astra_tcv_EC_exp(shot,[],[],[],[],[],gdat_data.gdat_params.trialindx); % centralized function for toray nodes + % centralized function for toray nodes + [pabs_gyro,icdtot,pow_dens,currentdrive_dens,rho_dep_pow,drho_pow,... + pmax,icdmax,currentdrive_dens_w2,rho_dep_icd,drho_icd,... + ~,power_integrated,currentdrive_integrated] = ... + astra_tcv_EC_exp(shot,[],[],[],[],[],gdat_data.gdat_params.trialindx,1); end if gdat_data.mapping_for.tcv.gdat_timedim ==2 tgrid_to_change = {'pabs_gyro','icdtot','pow_dens','currentdrive_dens','rho_dep_pow','drho_pow','pmax', ... @@ -1599,8 +1608,13 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end ec_help = 'from toray icdint with extracting of effective Icd for given launcher depending on nb rays used'; % All EC related quantities, each substructure should have at least fields data,x,t,units,dim,dimunits,label to be copied onto gdat_data - launchers_label = {'1','2','3','4','5','6','7','8','9','tot'}; - launchers_grid = [1:10]'; + if size(pabs_gyro.data,1)>9 + launchers_label = {'1','2','3','4','5','6','7','8','9','10','11','tot'}; + else + launchers_label = {'1','2','3','4','5','6','7','8','9','tot'}; + end + launchers_grid = [1:size(pabs_gyro.data,1)+1]'; + % power deposition related: ec_data.p_abs_plasma.data = pabs_gyro.data * 1e6; ec_data.p_abs_plasma.data(end+1,:) = sum(ec_data.p_abs_plasma.data,1,'omitnan'); @@ -1621,6 +1635,16 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ec_data.p_dens.dim = {ec_data.p_dens.x, launchers_grid, ec_data.p_dens.t}; ec_data.p_dens.dimunits = {'rhotor_norm', launchers_label, 's'}; % + ec_data.p_integrated.data = power_integrated.data * 1e6; + ec_data.p_integrated.data(:,end+1,:) = sum(ec_data.p_integrated.data,2,'omitnan'); + ec_data.p_integrated.label = [strrep(power_integrated.comment,'MW','W') ' ; last index is total']; + ec_data.p_integrated.units = 'W'; + ec_data.p_integrated.x = power_integrated.rgrid'; + ec_data.p_integrated.rhotor_norm = ec_data.p_integrated.x; + ec_data.p_integrated.t = power_integrated.tgrid; + ec_data.p_integrated.dim = {ec_data.p_integrated.x, launchers_grid, ec_data.p_integrated.t}; + ec_data.p_integrated.dimunits = {'rhotor_norm', launchers_label, 's'}; + % ec_data.max_pow_dens.data = pmax.data * 1e6; ec_data.max_pow_dens.label = strrep(pmax.comment,'MW','W'); ec_data.max_pow_dens.units = 'W/m^3'; @@ -1644,6 +1668,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ec_data.width_pow_dens.t = drho_pow.tgrid; ec_data.width_pow_dens.dim = {ec_data.width_pow_dens.t}; ec_data.width_pow_dens.dimunits = {'s'}; + % current drive deposition related: ec_data.cd_tot.data = icdtot.data * 1e6; ec_data.cd_tot.data(end+1,:) = sum(ec_data.cd_tot.data,1,'omitnan'); @@ -1662,7 +1687,17 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ec_data.cd_dens.rhotor_norm = ec_data.cd_dens.x; ec_data.cd_dens.t = currentdrive_dens.tgrid; ec_data.cd_dens.dim = {ec_data.cd_dens.x, launchers_grid, ec_data.cd_dens.t}; - ec_data.cd_dens.dimunits = {'rhotor_norm', launchers_label, 's'}; + ec_data.cd_dens.dimunits = {'rhotor_norm', launchers_label, 's'}; + % + ec_data.cd_integrated.data = currentdrive_integrated.data * 1e6; + ec_data.cd_integrated.data(:,end+1,:) = sum(ec_data.cd_integrated.data,2,'omitnan'); + ec_data.cd_integrated.label = [strrep(currentdrive_integrated.comment,'MA','A') ' ; last index is total']; + ec_data.cd_integrated.units = 'A'; + ec_data.cd_integrated.x = currentdrive_integrated.rgrid'; + ec_data.cd_integrated.rhotor_norm = ec_data.cd_integrated.x; + ec_data.cd_integrated.t = currentdrive_integrated.tgrid; + ec_data.cd_integrated.dim = {ec_data.cd_integrated.x, launchers_grid, ec_data.cd_integrated.t}; + ec_data.cd_integrated.dimunits = {'rhotor_norm', launchers_label, 's'}; % ec_data.max_cd_dens.data = icdmax.data * 1e6; ec_data.max_cd_dens.label = strrep(icdmax.comment,'MA','A'); @@ -1693,6 +1728,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') for subfields={'x','rhotor_norm','t','dim','dimunits','units'} ec_data.cd_dens_doublewidth.(subfields{1}) = ec_data.cd_dens.(subfields{1}); end + else disp(['source_icd.ec = ' source_icd.ec ' not yet implemented, ask O. Sauter']) ec_data.p_abs_plasma = []; @@ -1701,6 +1737,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ec_data.p_dens = []; ec_data.p_dens(end+1,:) = []; ec_data.p_dens_label = []; + ec_data.p_integrated = []; + ec_data.p_integrated(end+1,:) = []; + ec_data.p_integrated_label = []; ec_data.max_pow_dens = []; ec_data.max_pow_dens_label = []; ec_data.rho_max_pow_dens = []; @@ -1714,6 +1753,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ec_data.cd_dens = []; ec_data.cd_dens(:,end+1,:) = []; ec_data.cd_dens_label = []; + ec_data.cd_integrated = []; + ec_data.cd_integrated(:,end+1,:) = []; + ec_data.cd_integrated_label = []; ec_data.max_cd_dens = []; ec_data.max_cd_dens_label = []; ec_data.rho_max_cd_dens = []; @@ -3612,14 +3654,14 @@ catch gdat_data.error_bar_raw = gdat_data.error_bar; end - if (nverbose>=1) && shot<100000 + if (nverbose>=1) && shot<900000 warning('Problems with \results::thomson:times') warning(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff]) end return end if isempty(time) || ischar(time) - if (nverbose>=1) && shot<100000 + if (nverbose>=1) && shot<900000 warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check') disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff]) end diff --git a/matlab/TCV/loadTCVdata.m b/matlab/TCV/loadTCVdata.m deleted file mode 100644 index 49f83dcfe2d9e56ae0987a296c312b64e0808a23..0000000000000000000000000000000000000000 --- a/matlab/TCV/loadTCVdata.m +++ /dev/null @@ -1,1424 +0,0 @@ - function [trace,error_status,varargout]=loadTCVdata(shot,data_type,varargin) -% -% Added option to load shot=-1 or >=100000 -% Added option shot=-9 to list keywords -% -% list of data_type currently available (when [_2,_3] is added, means can add _i to get Liuqe i): -% if -1 is added, can also get it from FBTE with shot=-1, >=100000 or liuqe_version='_-1' (to get model file) -% -% 'Ip'[_2,_3] = current -% 'B0'[_2,_3] = current -% 'zmag'[_2,_3] = vertical position of the center of the plasma (magnetic axis) -% 'rmag'[_2,_3] = radial position of the center of the plasma -% 'rcont'[_2,_3] = R of plama boundary vs time -% 'zcont'[_2,_3] = Z of plama boundary vs time -% 'vol'[_2,_3] = volume of flux surfaces -% 'rhovol'[_2,_3] = sqrt(V(:,t)/V(edge,t)), normalised rho variable based on volume of flux surfaces -% 'qrho'[_2,_3] = q profile on rho mesh -% 'q95'[_2,_3] = q95 vs time -% 'kappa', 'elon'[_2,_3] = edge elongation vs time -% 'delta', 'triang'[_2,_3] = edge averaged triangularity vs time -% 'deltatop', 'triangtop'[_2,_3] = edge upper (top) triangularity vs time -% 'deltabot', 'triangbot'[_2,_3] = edge lower (bottom) triangularity vs time -% 'j_tor'[_2,_3] = J_TOR vs (R,Z,time) -% 'neint' = line-integrated electron density [m/m^3] -% 'nel' = line-averaged electron density [1/m^3] -% 'ne'= ne raw profile on (z,t). ADD error bars in .std -% 'te'= Te raw profile on (z,t). ADD error bars in .std -% 'nerho'[_2,_3]= ne profile on (rho=sqrt(psi),time) mesh.Note rho is a 2D array as depends on time. ADD error bars in .std -% 'terho'[_2,_3]= Te profile on (rho=sqrt(psi),time) mesh.Note rho is a 2D array as depends on time. ADD error bars in .std -% 'nerhozshift'[_2,_3]= same as nerho but allows for zshift [m] in equil given by varargin{1} -% 'terhozshift'[_2,_3]= same as terho but allows for zshift [m] in equil given by varargin{1} -% 'profnerho' = ne smoothed or fitted , vs (rho,t) (from Thomson auto fit) -% 'profterho' = te smoothed or fitted , vs (rho,t) (from Thomson auto fit) -% 'neft' = ne fitted from data on rho mesh (from proffit.local_time:neft_abs) -% 'neft:4' = ne fitted from data on rho mesh (from proffit.local_time:neft_abs from trial:trial_indx=4) -% 'teft' = te fitted from data on rho mesh (from proffit.local_time:teft) -% 'neftav' = ne fitted from averaged over time data on rho mesh (from proffit.avg_time:neft_abs) -% 'teftav' = te fitted from averaged over time data on rho mesh (from proffit.avg_time:teft) -% 'vtor'= Toroidal rotation of C raw profile on (rho,t). ADD error bars in .std, as well as time and rho error bars -% 'vtorfit'= Toroidal rotation of C fitted profile on (rho,t) -% 'vpol'= Poloidal rotation of C raw profile on (rho,t). ADD error bars in .std, as well as time and rho error bars -% 'vpolfit'= Poloidal rotation of C fitted profile on (rho,t) -% 'Ti' (or Tc), 'Tifit', 'ni' (or nC), 'nifit', 'zeffcxrs', 'zeffcxrsfit': similar to 'vtor' from CXRS -% 'ece' = electron cyclotron emission -% 'sxr' = soft x-ray emission -% 'sxR' = soft x-ray emission with varargout{1} option (requires varargin{4}!) -% 'MPX' = soft x-ray from wire chambers -% 'IOH' = current in ohmic transformer coils IOH_1 -% 'vloop' = loop voltage -% 'pgyro' = ECH power for each gyro(1:9) and total (10) -% -% 'xx_2 or xx_3' for Liuqe2 or 3: same as above for xx related to equilibrium -% 'xx_-1 or xx_-1' for FBTE values (model or shot=-1 or shot>=100000): same as above for xx related to equilibrium -% -% INPUT: -% shot: shot number -% data_type: type of the required data -% -% Definition of varargin depends on data_type: -% -% data_type=sxr or ece: -% varargin{1}: [i1 i2] : if not empty, assumes need many chords from i1 to i2 -% varargin{2}: channel status 1= unread yet, 0= read -% (for traces with many channel, enables to load additional channels, -% like SXR, ECE, etc.) -% varargin{3}: zmag for varargout{1} computation -% (can have more inputs for AUG, but not used here) -% data_type=nerhozshift or terhozshift: -% varargin{1}: zshift [m] constant or (t) : positive, moves equil up (that is thomson effective z down) -% time dependent: needs same time base as psitbx:psi -% -% OUTPUT: -% -% trace.data: data -% trace.t: time of reference -% trace.x: space of reference -% trace.dim: cell array of grids, trace.dim{1}, {2}, ... -% trace.dimunits: units of dimensions -% trace.units: units of data -% -% Additional Output arguments depending on data_type -% -% data_type=sxR, ece: -% varargout{1}: major radius: intersection/projection of the view lines with z=zmag -% data_type=MPX: -% varargout{1}: te be determined -% -% -% function needed: mds functions-xtomo_geometry-get_xtomo_data (Furno's routines) -% VsxrTCVradius -% Example: -% [ip,error_status]=loadTCVdata(shot,'Ip',1); -% [sxr,error_status,R]=loadTCVdata(shot,'sxR',1); -% - -varargout{1}=cell(1,1); -error_status=1; - -% To allow multiple ways of writing a specific keyword, use data_type_eff within this routine -if exist('data_type') && ~isempty(data_type) - data_type_eff=data_type; -else - data_type_eff=' '; -end -i_23=0; -% LIUQE tree -begstr = '\results::'; -endstr = ''; -liuqe_version = 1; -if length(data_type_eff)>2 - if strcmp(data_type_eff(end-1:end),'_2') | strcmp(data_type_eff(end-1:end),'_3') - i_23=2; - endstr=data_type_eff(end-1:end); - liuqe_version = str2num(data_type_eff(end:end)); - elseif strcmp(upper(data_type_eff(end-2:end)),'_-1') - i_23=3; - begstr = 'tcv_eq( "'; - endstr = '", "FBTE" )'; - liuqe_version = -1; - end -end -if shot==-1 || (shot>=100000 && shot<200000) - % requires FBTE - liuqe_version = -1; - begstr = 'tcv_eq( "'; - endstr = '", "FBTE" )'; -end - -% use keyword without eventual _2 or _3 extension to check for multiple possibilities -% Also remove ":4" for trial_indx specification -jj=strfind(data_type_eff,':'); -trialindx=[]; -if ~isempty(jj) - ii=strmatch(data_type_eff(1:jj-1),[{'teft'},{'neft'},{'teftav'},{'neftav'}]); - if ~isempty(ii) - trialindx=str2num(data_type_eff(jj+1:end)); - data_type_eff_noext=[data_type_eff(1:jj-1) ':trial']; - else - data_type_eff_noext=data_type_eff(1:end-i_23); - end -else - data_type_eff_noext=data_type_eff(1:end-i_23); -end - -if ~isempty(strmatch(data_type_eff_noext,[{'ip'} {'i_p'} {'xip'}],'exact')) - data_type_eff_noext='Ip'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'b0'} {'B_0'} {'Bgeom'} {'BGEOM'}],'exact')) - data_type_eff_noext='B0'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Te'} {'t_e'} {'TE'} {'T_e'}],'exact')) - data_type_eff_noext='te'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Te_edge'} {'TE_edge'}],'exact')) - data_type_eff_noext='te_edge'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Ne_edge'} {'NE_edge'}],'exact')) - data_type_eff_noext='ne'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Rcont'}],'exact')) - data_type_eff_noext='rcont'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Zcont'}],'exact')) - data_type_eff_noext='zcont'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Terho'}],'exact')) - data_type_eff_noext='terho'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Terhozshift'}],'exact')) - data_type_eff_noext='terhozshift'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'SXR'}],'exact')) - data_type_eff_noext='sxr'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'ECE'}],'exact')) - data_type_eff_noext='ece'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'VOL'} {'volume'}],'exact')) - data_type_eff_noext='vol'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'RHOVOL'}],'exact')) - data_type_eff_noext='rhovol'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'q_95'} {'Q95'}],'exact')) - data_type_eff_noext='q95'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'elongation'} {'elon'}],'exact')) - data_type_eff_noext='kappa'; -end -if ~isempty(strmatch(lower(data_type_eff_noext),[{'j_tor'} {'jtor'} {'\results::j_tor'}],'exact')) - data_type_eff_noext='jtor'; - data_type_eff = ['jtor' endstr]; -end -if ~isempty(strmatch(data_type_eff_noext,[{'triangularity'} {'triang'}],'exact')) - data_type_eff_noext='delta'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'deltaup'} {'deltau'} {'triangtop'} {'triangu'} {'triangup'}],'exact')) - data_type_eff_noext='deltatop'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'deltalow'} {'deltal'} {'triangbot'} {'triangl'} {'trianglow'}],'exact')) - data_type_eff_noext='deltabot'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Rmag'}],'exact')) - data_type_eff_noext='rmag'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Zmag'}],'exact')) - data_type_eff_noext='zmag'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'MPX'}],'exact')) - data_type_eff_noext='MPX'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'ioh'} {'Ioh'} {'iot'} {'IOT'}],'exact')) - data_type_eff_noext='IOH'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Vloop'} {'Vsurf'} {'vsurf'}],'exact')) - data_type_eff_noext='vloop'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'vtor'} {'v_tor'} {'vi'} {'vc'} {'v_i'} {'v_c'} {'VTOR'} {'V_TOR'} {'VI'} {'VC'} {'V_I'} {'V_C'}],'exact')) - data_type_eff_noext='vi_tor'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'vtorfit'} {'vtorft'} {'v_torfit'} {'vifit'} {'vcfit'} {'v_ifit'} {'v_cfit'} {'VTORfit'} {'V_TORfit'} {'VIfit'} {'VCfit'} {'V_Ifit'} {'V_Cfit'}],'exact')) - data_type_eff_noext='vi_torfit'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'vpol'} {'v_pol'} {'VPOL'} {'V_POL'}],'exact')) - data_type_eff_noext='vi_pol'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'vpolfit'} {'vpolft'} {'v_polfit'} {'VPOLfit'} {'V_POLfit'}],'exact')) - data_type_eff_noext='vi_polfit'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Ti'} {'ti'} {'TC'} {'tc'} {'T_C'} {'t_c'}],'exact')) - data_type_eff_noext='Ti'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Tifit'} {'tifit'} {'TCfit'} {'tcfit'} {'T_Cfit'} {'t_cfit'}],'exact')) - data_type_eff_noext='Tifit'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Ni'} {'ni'} {'NC'} {'nc'} {'N_C'} {'n_c'}],'exact')) - data_type_eff_noext='ni'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Nifit'} {'nifit'} {'NCfit'} {'ncfit'} {'N_Cfit'} {'n_cfit'}],'exact')) - data_type_eff_noext='nifit'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Zeffcxrs'} {'zeffcxrs'} {'Z_effcxrs'} {'z_effcxrs'} {'Zeff_cxrs'} {'zeff_cxrs'} {'Z_eff_cxrs'} {'z_eff_cxrs'}],'exact')) - data_type_eff_noext='zeffcxrs'; -end -if ~isempty(strmatch(data_type_eff_noext,[{'Zeffcxrsfit'} {'zeffcxrsfit'} {'Z_effcxrsfit'} {'z_effcxrsfit'} {'Zeff_cxrsfit'} {'zeff_cxrsfit'} {'Z_eff_cxrsfit'} {'z_eff_cxrsfit'}],'exact')) - data_type_eff_noext='zeffcxrsfit'; -end - -% some defaults - -% nodes which have _2 and _3 equivalence, related to simpletdi case -liuqe23=[{'\results::area'} {'\results::beta_pol'} {'\results::beta_tor'} ... - {'\results::delta_95'} {'\results::delta_95_bot'} {'\results::delta_95_top'} ... - {'\results::delta_edge'} {'\results::delta_ed_bot'} {'\results::delta_ed_top'} ... - {'\results::diamag'} {'\results::i_p'} {'\results::j_tor'} ... - {'\results::kappa_95'} {'\results::kappa_edge'} {'\results::kappa_zero'} ... - {'\results::lambda'} {'\results::l_i'} {'\results::psi_axis'} {'\results::psi_values'} ... - {'\results::q_95'} {'\results::q_edge'} {'\results::q_zero'} {'\results::rms_error'} ... - {'\results::z_axis'} {'\results::r_axis'} {'\results::q_psi'} ... - {'\results::total_energy'} {'\results::tor_flux'} {'\results::volume'} ... - {'\results::r_contour'} {'\results::z_contour'} ... - {'\results::thomson:psiscatvol'} {'\results::thomson:psi_max'}]; - -% nodes which have FBTE equivalence, related to simpletdi case -liuqeFBTE=[{'\results::i_p'} {'\results::z_axis'} {'\results::r_axis'} {'\results::q_psi'} ... - {'\results::beta_tor'} {'\results::beta_pol'} {'\results::q_95'} {'\results::l_i'} {'\results::delta_95'} ... - {'\results::kappa_95'} {'\results::r_contour'} {'\results::z_contour'} {'\results::psi_axis'}]; - -% keywords which do not have FBTE equivalence and are returned empty -noFBTE=[{'ne'} {'te'} {'nerho'} {'terho'} {'ne_edge'} {'te_edge'} {'nerho_edge'} {'terho_edge'} {'nerhozshift'} {'terhozshift'} {'profnerho'} {'profterho'} ... - {'neft'} {'neft:trial'} {'teft:trial'} {'neftav:trial'} {'teftav:trial'} {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'} {'neint'} {'nel'} ... - {'vi_tor'} {'vi_torfit'} {'vi_pol'} {'vi_polfit'} {'Ti'} {'Tifit'} {'ni'} {'nifit'} {'zeffcxrs'} {'zeffcxrsfit'}]; - -% all keywords and corresponding case to run below -TCVkeywrdall=[{'Ip'} {'B0'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'rhovol'} {'qrho'} {'q95'} {'kappa'} ... - {'delta'} {'deltatop'} {'deltabot'} {'neint'} {'nel'} ... - {'ne'} {'te'} {'nerho'} {'terho'} {'ne_edge'} {'te_edge'} {'nerho_edge'} {'terho_edge'} {'nerhozshift'} {'terhozshift'} {'profnerho'} {'profterho'} ... - {'neft'} {'teft'} {'neftav'} {'teftav'} {'neft:trial'} {'teft:trial'} {'neftav:trial'} {'teftav:trial'} ... - {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'} {'pgyro'} {'jtor'} {'vi_tor'} {'vi_torfit'} {'vi_pol'} {'vi_polfit'} {'Ti'} {'Tifit'} {'ni'} {'nifit'} {'zeffcxrs'} {'zeffcxrsfit'}]; - -TCVsig.iip=strmatch('Ip',TCVkeywrdall,'exact'); -TCVsig.iB0=strmatch('B0',TCVkeywrdall,'exact'); -TCVsig.izmag=strmatch('zmag',TCVkeywrdall,'exact'); -TCVsig.irmag=strmatch('rmag',TCVkeywrdall,'exact'); -TCVsig.ircont=strmatch('rcont',TCVkeywrdall,'exact'); -TCVsig.izcont=strmatch('zcont',TCVkeywrdall,'exact'); -TCVsig.ivol=strmatch('vol',TCVkeywrdall,'exact'); -TCVsig.irhovol=strmatch('rhovol',TCVkeywrdall,'exact'); -TCVsig.iqrho=strmatch('qrho',TCVkeywrdall,'exact'); -TCVsig.iq95=strmatch('q95',TCVkeywrdall,'exact'); -TCVsig.ikappa=strmatch('kappa',TCVkeywrdall,'exact'); -TCVsig.idelta=strmatch('delta',TCVkeywrdall,'exact'); -TCVsig.ideltatop=strmatch('deltatop',TCVkeywrdall,'exact'); -TCVsig.ideltabot=strmatch('deltabot',TCVkeywrdall,'exact'); -TCVsig.ineint=strmatch('neint',TCVkeywrdall,'exact'); -TCVsig.inel=strmatch('nel',TCVkeywrdall,'exact'); -TCVsig.ine=strmatch('ne',TCVkeywrdall,'exact'); -TCVsig.ite=strmatch('te',TCVkeywrdall,'exact'); -TCVsig.ine_edge=strmatch('ne_edge',TCVkeywrdall,'exact'); -TCVsig.ite_edge=strmatch('te_edge',TCVkeywrdall,'exact'); -TCVsig.inerho=strmatch('nerho',TCVkeywrdall,'exact'); -TCVsig.iterho=strmatch('terho',TCVkeywrdall,'exact'); -TCVsig.inerho_edge=strmatch('nerho_edge',TCVkeywrdall,'exact'); -TCVsig.iterho_edge=strmatch('terho_edge',TCVkeywrdall,'exact'); -TCVsig.inerhozshift=strmatch('nerhozshift',TCVkeywrdall,'exact'); -TCVsig.iterhozshift=strmatch('terhozshift',TCVkeywrdall,'exact'); -TCVsig.iprofnerho=strmatch('profnerho',TCVkeywrdall,'exact'); -TCVsig.iprofterho=strmatch('profterho',TCVkeywrdall,'exact'); -TCVsig.ineft=strmatch('neft',TCVkeywrdall,'exact'); -TCVsig.ineft_trial=strmatch('neft:trial',TCVkeywrdall); -TCVsig.iteft=strmatch('teft',TCVkeywrdall,'exact'); -TCVsig.iteft_trial=strmatch('teft:trial',TCVkeywrdall); -TCVsig.ineftav=strmatch('neftav',TCVkeywrdall,'exact'); -TCVsig.ineftav_trial=strmatch('neftav:trial',TCVkeywrdall); -TCVsig.iteftav=strmatch('teftav',TCVkeywrdall,'exact'); -TCVsig.iteftav_trial=strmatch('teftav:trial',TCVkeywrdall); -TCVsig.isxr=strmatch('sxr',TCVkeywrdall,'exact'); -TCVsig.isxR=strmatch('sxR',TCVkeywrdall,'exact'); -TCVsig.iece=strmatch('ece',TCVkeywrdall,'exact'); -TCVsig.iMPX=strmatch('MPX',TCVkeywrdall,'exact'); -TCVsig.iIOH=strmatch('IOH',TCVkeywrdall,'exact'); -TCVsig.ivloop=strmatch('vloop',TCVkeywrdall,'exact'); -TCVsig.ipgyro=strmatch('pgyro',TCVkeywrdall,'exact'); -TCVsig.ijtor=strmatch('jtor',TCVkeywrdall,'exact'); -TCVsig.ivi_tor=strmatch('vi_tor',TCVkeywrdall,'exact'); -TCVsig.ivi_torfit=strmatch('vi_torfit',TCVkeywrdall,'exact'); -TCVsig.ivi_pol=strmatch('vi_pol',TCVkeywrdall,'exact'); -TCVsig.ivi_polfit=strmatch('vi_polfit',TCVkeywrdall,'exact'); -TCVsig.iTi=strmatch('Ti',TCVkeywrdall,'exact'); -TCVsig.iTifit=strmatch('Tifit',TCVkeywrdall,'exact'); -TCVsig.ini=strmatch('ni',TCVkeywrdall,'exact'); -TCVsig.inifit=strmatch('nifit',TCVkeywrdall,'exact'); -TCVsig.izeffcxrs=strmatch('zeffcxrs',TCVkeywrdall,'exact'); -TCVsig.izeffcxrsfit=strmatch('zeffcxrsfit',TCVkeywrdall,'exact'); - -% For each keyword, specify which case to use. As most common is 'simpletdi', fill in with this and change -% only indices needed. Usually use name of case same as keyword name -TCVkeywrdcase=cell(size(TCVkeywrdall)); -TCVkeywrdcase(:)={'simpletdi'}; -TCVkeywrdcase(TCVsig.iB0)=TCVkeywrdall(TCVsig.iB0); % through iphi -TCVkeywrdcase(TCVsig.iqrho)=TCVkeywrdall(TCVsig.iqrho); % special as liuqe q_psi on psi -TCVkeywrdcase(TCVsig.ivol)=TCVkeywrdall(TCVsig.ivol); % special as nodes _2 or _3 not existing with psitbx -TCVkeywrdcase(TCVsig.irhovol)=TCVkeywrdall(TCVsig.irhovol); % idem vol -TCVkeywrdcase(TCVsig.ine)=TCVkeywrdall(TCVsig.ine); % special as dimensions from other nodes -TCVkeywrdcase(TCVsig.ite)=TCVkeywrdall(TCVsig.ite); % idem -TCVkeywrdcase(TCVsig.ine_edge)=TCVkeywrdall(TCVsig.ine_edge); % special as dimensions from other nodes -TCVkeywrdcase(TCVsig.ite_edge)=TCVkeywrdall(TCVsig.ite_edge); % idem -TCVkeywrdcase(TCVsig.inerho)=TCVkeywrdall(TCVsig.inerho); % idem -TCVkeywrdcase(TCVsig.inerho_edge)=TCVkeywrdall(TCVsig.inerho_edge); % idem -TCVkeywrdcase(TCVsig.iteft_trial)=TCVkeywrdall(TCVsig.iteft_trial); % special to extract trial_indx if needed -TCVkeywrdcase(TCVsig.ineft_trial)=TCVkeywrdall(TCVsig.ineft_trial); % special to extract trial_indx if needed -TCVkeywrdcase(TCVsig.iteftav_trial)=TCVkeywrdall(TCVsig.iteftav_trial); % special to extract trial_indx if needed -TCVkeywrdcase(TCVsig.ineftav_trial)=TCVkeywrdall(TCVsig.ineftav_trial); % special to extract trial_indx if needed -TCVkeywrdcase(TCVsig.iterho)=TCVkeywrdall(TCVsig.iterho); % idem -TCVkeywrdcase(TCVsig.iterho_edge)=TCVkeywrdall(TCVsig.iterho_edge); % idem -TCVkeywrdcase(TCVsig.inerhozshift)=TCVkeywrdall(TCVsig.inerhozshift); % idem -TCVkeywrdcase(TCVsig.iterhozshift)=TCVkeywrdall(TCVsig.iterhozshift); % idem -TCVkeywrdcase(TCVsig.iprofnerho)=TCVkeywrdall(TCVsig.iprofnerho); -TCVkeywrdcase(TCVsig.iprofterho)=TCVkeywrdall(TCVsig.iprofterho); -TCVkeywrdcase(TCVsig.isxr)=TCVkeywrdall(TCVsig.isxr); -TCVkeywrdcase(TCVsig.isxR)=TCVkeywrdall(TCVsig.isxR); -TCVkeywrdcase(TCVsig.iece)=TCVkeywrdall(TCVsig.iece); -TCVkeywrdcase(TCVsig.iMPX)=TCVkeywrdall(TCVsig.iMPX); -TCVkeywrdcase(TCVsig.iIOH)=TCVkeywrdall(TCVsig.iIOH); -TCVkeywrdcase(TCVsig.ivloop)=TCVkeywrdall(TCVsig.ivloop); -TCVkeywrdcase(TCVsig.ipgyro)=TCVkeywrdall(TCVsig.ipgyro); -TCVkeywrdcase(TCVsig.ijtor)=TCVkeywrdall(TCVsig.ijtor); -TCVkeywrdcase(TCVsig.ivi_tor)=TCVkeywrdall(TCVsig.ivi_tor); -TCVkeywrdcase(TCVsig.ivi_torfit)=TCVkeywrdall(TCVsig.ivi_torfit); -TCVkeywrdcase(TCVsig.ivi_pol)=TCVkeywrdall(TCVsig.ivi_pol); -TCVkeywrdcase(TCVsig.ivi_polfit)=TCVkeywrdall(TCVsig.ivi_polfit); -TCVkeywrdcase(TCVsig.iTi)=TCVkeywrdall(TCVsig.iTi); -TCVkeywrdcase(TCVsig.iTifit)=TCVkeywrdall(TCVsig.iTifit); -TCVkeywrdcase(TCVsig.ini)=TCVkeywrdall(TCVsig.ini); -TCVkeywrdcase(TCVsig.inifit)=TCVkeywrdall(TCVsig.inifit); -TCVkeywrdcase(TCVsig.izeffcxrs)=TCVkeywrdall(TCVsig.izeffcxrs); -TCVkeywrdcase(TCVsig.izeffcxrsfit)=TCVkeywrdall(TCVsig.izeffcxrsfit); - -% Information about which dimension has time, always return 2D data as (x,t) array -% as most are 1D arrays with time as first index, fill in with ones and change only those needed -TCVsigtimeindx=ones(size(TCVkeywrdall)); - -% For the 'simpletdi' cases, we need the full node in case gdat was called with full location directly -% for the other cases, leave this location empty -TCVsiglocation=cell(size(TCVkeywrdall)); -TCVsiglocation(:)={''}; -TCVsiglocation(TCVsig.iip)={'\results::i_p'}; -TCVsiglocation(TCVsig.izmag)={'\results::z_axis'}; -TCVsiglocation(TCVsig.irmag)={'\results::r_axis'}; -TCVsiglocation(TCVsig.ircont)={'\results::r_contour'}; TCVsigtimeindx(TCVsig.ircont)=2; -TCVsiglocation(TCVsig.izcont)={'\results::z_contour'}; TCVsigtimeindx(TCVsig.izcont)=2; -TCVsiglocation(TCVsig.iq95)={'\results::q_95'}; -TCVsiglocation(TCVsig.ikappa)={'\results::kappa_edge'}; -TCVsiglocation(TCVsig.idelta)={'\results::delta_edge'}; -TCVsiglocation(TCVsig.ideltatop)={'\results::delta_ed_top'}; -TCVsiglocation(TCVsig.ideltabot)={'\results::delta_ed_bot'}; -TCVsiglocation(TCVsig.ineint)={'\results::fir:lin_int_dens'}; -TCVsiglocation(TCVsig.inel)={'\results::fir:n_average'}; -%TCVsiglocation(TCVsig.iprofnerho)={'\results::THOMSON.PROFILES.AUTO:ne'}; -%TCVsiglocation(TCVsig.iprofterho)={'\results::THOMSON.PROFILES.AUTO:te'}; -TCVsiglocation(TCVsig.ineft)={'\results::proffit.local_time:neft_abs'}; TCVsigtimeindx(TCVsig.ineft)=2; -TCVsiglocation(TCVsig.iteft)={'\results::proffit.local_time:teft'}; TCVsigtimeindx(TCVsig.iteft)=2; -TCVsiglocation(TCVsig.ineftav)={'\results::proffit.avg_time:neft_abs'}; TCVsigtimeindx(TCVsig.ineftav)=2; -TCVsiglocation(TCVsig.iteftav)={'\results::proffit.avg_time:teft'}; TCVsigtimeindx(TCVsig.iteftav)=2; -TCVsiglocation(TCVsig.ineft_trial)=TCVsiglocation(TCVsig.ineft); TCVsigtimeindx(TCVsig.ineft_trial)=2; -TCVsiglocation(TCVsig.iteft_trial)=TCVsiglocation(TCVsig.iteft); TCVsigtimeindx(TCVsig.iteft_trial)=2; -TCVsiglocation(TCVsig.ineftav_trial)=TCVsiglocation(TCVsig.ineftav); TCVsigtimeindx(TCVsig.ineftav_trial)=2; -TCVsiglocation(TCVsig.iteftav_trial)=TCVsiglocation(TCVsig.iteftav); TCVsigtimeindx(TCVsig.iteftav_trial)=2; - -if shot==-9 - clear trace - for i=1:length(TCVkeywrdall) - fieldname_eff = lower(TCVkeywrdall{i}); - keyword_eff = TCVkeywrdall{i}; - ij=findstr(fieldname_eff,':'); - if ~isempty(ij); - fieldname_eff(ij)='_'; - keyword_eff(ij:end+1)=[':i' keyword_eff(ij+1:end)] ; - end - trace.(fieldname_eff) = keyword_eff; - end - % add example for Ip trace full call - trace.ipfullcall = TCVsiglocation{TCVsig.iip}; - % trace.data=[]; - return -end - -% initialize order of substructures and allows just a "return" if data empty -trace.data=[]; -trace.x=[]; -trace.t=[]; -trace.dim=[]; -trace.dimunits=[]; -trace.units=[]; -iprintwarn=0; -% find index of signal called upon -if strcmp(data_type_eff(1:1),'\') - % in case full node name was given - index=strmatch(data_type_eff(1:end-i_23),TCVsiglocation,'exact'); - if iprintwarn & isempty(index) - disp('********************') - disp('trace not yet registered.') - disp('If standard data, ask olivier.sauter@epfl.ch to create a keyqord entry for this data') -% eval(['!mail -s ''' data_type_eff ' ' num2str(shot) ' ' getenv('USER') ' TCV'' olivier.sauter@epfl.ch < /dev/null']) - disp('********************') - elseif isempty(index) - % temporarily add entry in arrays, so can work below - index=length(TCVkeywrdall)+1; - TCVkeywrdall(end+1)={'new'}; - TCVkeywrdcase(end+1)={'simpletdi'}; - TCVsiglocation(end+1)={data_type_eff(1:end-i_23)}; - TCVsigtimeindx(end+1)=0; - elseif ~strcmp(TCVkeywrdcase{index},'simpletdi') - msgbox(['Problem in loadTCVdata with data_type_eff = ' data_type_eff ... - '. Full paths of nodes should only be for case simpletdi'],'in loadTCVdata','error') - error('in loadTCVdata') - end -else - index=strmatch(data_type_eff_noext,TCVkeywrdall,'exact'); - if isempty(index) - disp(' ') - disp('********************') - disp(['no such keyword: ' data_type_eff(1:end-i_23)]) - disp(' ') - disp('Available keywords:') - TCVkeywrdall(:) - disp('********************') - return - end - if liuqe_version==-1 && ~isempty(strmatch(TCVkeywrdall{index},noFBTE,'exact')) - disp(['node ' TCVkeywrdall{index} ' not defined within FBTE']) - return - end -end -if iprintwarn - disp(['loading' ' ' data_type_eff_noext ' from TCV shot #' num2str(shot)]); - disp(['case ' TCVkeywrdcase{index}]) -end -if i_23==2 & isempty(strmatch(TCVsiglocation(index),liuqe23,'exact')) & strcmp(TCVkeywrdcase{index},'simpletdi') - disp('********') - disp('Warning asks for liuqe 2 or 3 of a signal, but not in liuqe23 list in loadTCVdata') -% eval(['!mail -s ''' data_type_eff ' ' num2str(shot) ' ' getenv('USER') ' TCV: not in liuqe23 list'' olivier.sauter@epfl.ch < /dev/null']) - disp('********') -end - -status=ones(1,100); -zmag=cell(0,0); -nargineff=nargin; -if iprintwarn - disp('TCVsiglocation{index} in loadTCVdata') - TCVsiglocation{index} -end - -switch TCVkeywrdcase{index} - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case 'simpletdi' - - % load TCV other data - if liuqe_version==-1 - mdsopen( 'pcs', shot); %synthetic shot generated with FBTE and MGAMS. - else - mdsopen(shot); - % test if node exists - error_status=1; - ij=findstr(TCVsiglocation{index},'['); - if isempty(ij); ij=length(TCVsiglocation{index})+1; end - if eval(['~mdsdata(''node_exists("\' TCVsiglocation{index}(1:ij-1) '")'')']) - disp(['node ' TCVsiglocation{index}(1:ij-1) ' does not exist for shot = ' num2str(shot)]) - return -% $$$ elseif eval(['mdsdata(''getnci("\' nodenameeff ':foo","length")'')==0']) -% $$$ disp(['no data for node ' nodenameeff ' for shot = ' num2str(shot)]) -% $$$ return - end - end - if strcmp(TCVsiglocation{index}(1:9),'\psitbx::') - begstr = 'tcv_psitbx("'; - nodenameeff=[begstr TCVsiglocation{index}(10:end) endstr]; - elseif length(TCVsiglocation{index})>16 && strcmp(TCVsiglocation{index}(1:16),'\results::psitbx') - begstr = 'tcv_psitbx("'; - nodenameeff=[begstr TCVsiglocation{index}(18:end) endstr]; - elseif strcmp(TCVsiglocation{index}(1:10),'\results::') - nodenameeff=[begstr TCVsiglocation{index}(11:end) endstr]; - else - nodenameeff=TCVsiglocation{index}; - disp(['should not have gone through here, mention this example to O. Sauter']); - end - tracetdi=tdi(nodenameeff); - mdsclose; - if isempty(tracetdi.data) || (~iscell(tracetdi.data) && (length(tracetdi.data)<=1) && isnan(tracetdi.data)) - disp(['node ' nodenameeff ' is empty for shot = ' num2str(shot)]) - return - end - trace.data=tracetdi.data; - - if length(tracetdi.dim)==0 - % scalar - trace.x=[]; - trace.t=[]; - elseif TCVsigtimeindx(index)>0 | length(tracetdi.dim)==1 - trace.t=tracetdi.dim{max(1,TCVsigtimeindx(index))}; - ix=3-TCVsigtimeindx(index); % works only for 2D arrays - else - % if dim=41 assumes is x (usual for TCV) - if length(tracetdi.dim{1})==41 - ix=1; - trace.t=tracetdi.dim{2}; - elseif length(tracetdi.dim{2})==41 - ix=2; - trace.t=tracetdi.dim{1}; - elseif length(tracetdi.dim)==2 - % if one has 2D or more and the other 1D, assume rho is 2D or more and time is 1D - if min(size(tracetdi.dim{1}))==1 && (min(size(tracetdi.dim{2}))>1 || length(size(tracetdi.dim{2}))>2) - ix=2; - trace.t=tracetdi.dim{1}; - disp(['assumes time array is 1D array, so x from dim{' num2str(ix) '}']) - elseif min(size(tracetdi.dim{2}))==1 && (min(size(tracetdi.dim{1}))>1 || length(size(tracetdi.dim{1}))>2) - ix=1; - trace.t=tracetdi.dim{2}; - disp(['assumes time array is 1D array, so x from dim{' num2str(ix) '}']) - elseif tracetdi.dim{1}(end)<1.1 && tracetdi.dim{2}(end)>1.1 - ix=1; - trace.t=tracetdi.dim{2}; - disp(['assumes x array is 1D array with max <1.1']) - elseif tracetdi.dim{1}(end)>1.1 && tracetdi.dim{2}(end)<1.1 - ix=2; - trace.t=tracetdi.dim{1}; - disp(['assumes x array is 1D array with max <1.1']) - elseif length(size(tracetdi.dim{1}))==2 && min(size(tracetdi.dim{1}))==1 && length(size(tracetdi.dim{2}))==2 && min(size(tracetdi.dim{2}))==1 - % if both are 1D arrays, assumes time array with values having most number of digits - ab1=num2str(tracetdi.dim{1}-fix(tracetdi.dim{1})); - ab2=num2str(tracetdi.dim{2}-fix(tracetdi.dim{2})); - if size(ab1,2)<size(ab2,2); - ix=1; - trace.t=tracetdi.dim{2}; - else - ix=2; - trace.t=tracetdi.dim{1}; - end - disp(['assumes time array has more digits, so x from dim{' num2str(ix) '}']) - end - end - end - if length(tracetdi.dim)==2 - trace.x=tracetdi.dim{ix}; - % make sure data is of (x,t) - if ix==2 - % transpose data - trace.data=trace.data'; - end - elseif length(tracetdi.dim)>2 - msgbox(['data more than 2D (data_type_eff=' data_type_eff ... - '), check how to save it, contact andrea.scarabosio@epfl.ch or olivier.sauter@epfl.ch'],... - 'in simpletdi','warn') - warning('in simpletdi of loadTCVdata') - else - if max(1,TCVsigtimeindx(index))~=1 - disp('Problems in loadTCVdata, max(1,TCVsigtimeindx(index)) should be 1') - end - end - if length(tracetdi.dim)==1 - trace.dim=[{trace.t}]; - trace.dimunits=tracetdi.dimunits(1); - elseif length(tracetdi.dim)==2 - trace.dim=[{trace.x} ; {trace.t}]; - trace.dimunits=[tracetdi.dimunits(ix) ; tracetdi.dimunits((2-ix)+1)]; - else - trace.dim=tracetdi.dim; - trace.dimunits=tracetdi.dimunits; - end - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - - trace.name=[num2str(shot) ';' nodenameeff]; - error_status=0; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'ne','te'} - % ne or Te from Thomson data on raw z mesh vs (z,t) - mdsopen(shot); - if strcmp(TCVkeywrdcase{index},'ne') - nodenameeff='\results::thomson:ne'; - tracetdi=tdi(nodenameeff); - tracestd=tdi('\results::thomson:ne:error_bar'); - else - nodenameeff='\results::thomson:te'; - tracetdi=tdi(nodenameeff); - tracestd=tdi('\results::thomson:te:error_bar'); - trace_abs=tdi('\results::thomson:fir_thom_rat'); - end - trace.data=tracetdi.data'; % Thomson data as (t,z) - trace.std=tracestd.data'; - trace.name=[num2str(shot) ';' nodenameeff]; - % add correct dimensions - try - time=mdsdata('\results::thomson:times'); - catch - warning('Problems with \results::thomson:times') - disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}]) - return - end - if isempty(time) || ischar(time) - thomsontimes=time - warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check') - disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}]) - return - end - z=mdsdata('\diagz::thomson_set_up:vertical_pos'); - trace.dim=[{z};{time}]; - trace.dimunits=[{'Z [m]'} ; {'time [s]'}]; - trace.x=z; - trace.t=time; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - mdsclose('mdsip'); - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'ne_edge','te_edge'} - % ne or Te from Thomson.edge data on raw z mesh vs (z,t) - mdsopen(shot); - if strcmp(TCVkeywrdcase{index},'ne_edge') - nodenameeff='\results::thomson.edge:ne'; - tracetdi=tdi(nodenameeff); - tracestd=tdi('\results::thomson.edge:ne:error_bar'); - else - nodenameeff='\results::thomson.edge:te'; - tracetdi=tdi(nodenameeff); - tracestd=tdi('\results::thomson.edge:te:error_bar'); - trace_abs=tdi('\results::thomson.edge:fir_thom_rat'); - end - trace.data=tracetdi.data'; % Thomson.Edge data as (t,z) - trace.std=tracestd.data'; - trace.name=[num2str(shot) ';' nodenameeff]; - % add correct dimensions - try - time=mdsdata('\results::thomson:times'); - catch - warning('Problems with \results::thomson:times') - disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}]) - return - end - if isempty(time) || ischar(time) - thomsontimes=time - warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check') - disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}]) - return - end - z=mdsdata('\diagz::thomson_set_up.edge:vertical_pos'); - trace.dim=[{z};{time}]; - trace.dimunits=[{'Z [m]'} ; {'time [s]'}]; - trace.x=z; - trace.t=time; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - mdsclose('mdsip'); - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'nerho','terho'} - % ne or Te from Thomson data on rho=sqrt(psi_normalised) mesh: (rho,t) - mdsopen(shot); - try - time=mdsdata('\results::thomson:times'); - catch - warning('Problems with \results::thomson:times') - disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}]) - return - end - if isempty(time) || ischar(time) - thomsontimes=time - warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check') - disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}]) - return - end - if strcmp(TCVkeywrdcase{index},'nerho') - nodenameeff='\results::thomson:ne'; - tracetdi=tdi(nodenameeff); - if isempty(tracetdi.data) - ishot=mdsopen(shot) - tracetdi=tdi(nodenameeff) - end - nodenameeff='\results::thomson:ne; error_bar ; fir_thom_rat; (ne,std)*fir_thom_rat'; - tracestd=tdi('\results::thomson:ne:error_bar'); - if shot>=23801 - tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!! - if isempty(tracefirrat.data) || ischar(tracefirrat.data) - disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty') - end - else - tracefirrat=tdi('\results::thomson:fir_thom_rat'); - if isempty(tracefirrat.data) || ischar(tracefirrat.data) - disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty') - tracefirrat.dim{1}=[]; - else - tracefirrat.dim{1}=time; - end - end - if ~isempty(tracefirrat.data) || ischar(tracefirrat.data) - tracefirrat_data=NaN*ones(size(tracetdi.dim{1})); - itim=iround(time,tracefirrat.dim{1}); - tracefirrat_data(itim)=tracefirrat.data; - else - tracefirrat_data=NaN; - end - else - nodenameeff='\results::thomson:te'; - tracetdi=tdi(nodenameeff); - nodenameeff='\results::thomson:te; error_bar'; - tracestd=tdi('\results::thomson:te:error_bar'); - end - trace.data=tracetdi.data'; % Thomson data as (t,z) - trace.std=tracestd.data'; - if strcmp(TCVkeywrdcase{index},'nerho') - trace.firrat=tracefirrat_data; - trace.data_abs=trace.data*diag(tracefirrat_data); - trace.std_abs=trace.std*diag(tracefirrat_data); - end - trace.name=[num2str(shot) ';' nodenameeff]; - % add correct dimensions - % construct rho mesh - psi_max=tdi(['\results::thomson:psi_max' endstr]); - psiscatvol=tdi(['\results::thomson:psiscatvol' endstr]); - if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data) - for ir=1:length(psiscatvol.dim{2}) - rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))'; - end - else - rho=NaN; - end - trace.dim=[{rho};{time}]; - trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}]; - trace.x=rho; - trace.t=time; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - mdsclose; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'nerho_edge','terho_edge'} - % ne or Te from Thomson.Edge data on rho=sqrt(psi_normalised) mesh: (rho,t) - mdsopen(shot); - try - time=mdsdata('\results::thomson:times'); - catch - warning('Problems with \results::thomson:times') - disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}]) - return - end - if isempty(time) || ischar(time) - thomsontimes=time - warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check') - disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}]) - return - end - if strcmp(TCVkeywrdcase{index},'nerho_edge') - nodenameeff='\results::thomson.edge:ne'; - tracetdi=tdi(nodenameeff); - if isempty(tracetdi.data) || ischar(tracetdi.data) - ishot=mdsopen(shot) - tracetdi=tdi(nodenameeff) - end - nodenameeff='\results::thomson.edge:ne; error_bar ; fir_thom_rat; (ne,std)*fir_thom_rat'; - tracestd=tdi('\results::thomson.edge:ne:error_bar'); - if shot>=23801 - tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!! - if isempty(tracefirrat.data) || ischar(tracefirrat.data) - disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty') - end - else - tracefirrat=tdi('\results::thomson.edge:fir_thom_rat'); - tracefirrat.dim{1}=time; - end - if ~isempty(tracetdi.dim) - tracefirrat_data=NaN*ones(size(tracetdi.dim{1})); - else - tracefirrat_data=[]; - end - if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) - itim=iround(time,tracefirrat.dim{1}); - tracefirrat_data(itim)=tracefirrat.data; - end - else - nodenameeff='\results::thomson.edge:te'; - tracetdi=tdi(nodenameeff); - nodenameeff='\results::thomson.edge:te; error_bar'; - tracestd=tdi('\results::thomson.edge:te:error_bar'); - end - if ~ischar(tracetdi.data); trace.data=tracetdi.data'; end % Thomson.Edge data as (t,z) - trace.std=tracestd.data'; - if strcmp(TCVkeywrdcase{index},'nerho') - trace.firrat=tracefirrat_data; - trace.data_abs=trace.data*diag(tracefirrat_data); - trace.std_abs=trace.std*diag(tracefirrat_data); - end - trace.name=[num2str(shot) ';' nodenameeff]; - % add correct dimensions - % construct rho mesh - psi_max=tdi(['\results::thomson.edge:psi_max' endstr]); - psiscatvol=tdi(['\results::thomson.edge:psiscatvol' endstr]); - if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) - for ir=1:length(psiscatvol.dim{2}) - rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))'; - end - trace.dim=[{rho};{time}]; - trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}]; - trace.x=rho; - trace.t=time; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - else - trace.dim={}; - trace.dimunits={}; - trace.x=[]; - trace.t=[]; - trace.units={}; - end - mdsclose; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'nerhozshift','terhozshift'} - % ne or Te from Thomson data on rho=sqrt(psi_normalised) mesh: (rho,t) - % allow for z shift of equil - if nargin>=3 & ~isempty(varargin{1}) - zshift=varargin{1}; - else - zshift=0.; - end - mdsopen(shot); - try - time=mdsdata('\results::thomson:times'); - catch - warning('Problems with \results::thomson:times') - disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}]) - return - end - if isempty(time) || ischar(time) - thomsontimes=time - warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times is empty? Check') - disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}]) - return - end - if strcmp(TCVkeywrdcase{index},'nerhozshift') - nodenameeff='\results::thomson:ne'; - tracetdi=tdi(nodenameeff); - tracestd=tdi('\results::thomson:ne:error_bar'); - if shot>=23801 - tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!! - if isempty(tracefirrat.data) || ischar(tracefirrat.data) - disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty') - end - else - tracefirrat=tdi('\results::thomson:fir_thom_rat'); - if isempty(tracefirrat.data) || ischar(tracefirrat.data) - disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty') - tracefirrat.dim{1}=[]; - else - tracefirrat.dim{1}=time; - end - end - if ~isempty(tracefirrat.data) || ischar(tracefirrat.data) - tracefirrat_data=NaN*ones(size(tracetdi.dim{1})); - itim=iround(time,tracefirrat.dim{1}); - tracefirrat_data(itim)=tracefirrat.data; - else - tracefirrat_data=NaN; - end - else - nodenameeff='\results::thomson:te'; - tracetdi=tdi(nodenameeff); - tracestd=tdi('\results::thomson:te:error_bar'); - end - trace.data=tracetdi.data'; % Thomson data as (t,z) - trace.std=tracestd.data'; - if strcmp(TCVkeywrdcase{index},'nerhozshift') - trace.firrat=tracefirrat_data; - trace.data_abs=trace.data*diag(tracefirrat_data); - trace.std_abs=trace.std*diag(tracefirrat_data); - end - trace.name=[num2str(shot) ';' nodenameeff]; - % add correct dimensions - % construct rho mesh - if strcmp(endstr,'_-1') - error(['in ' TCVkeywrdcase{index} ' endstr should not be ' endstr]); - end - psi_max=tdi(['\results::thomson:psi_max' endstr]); - psiscatvol=tdi(['\results::thomson:psiscatvol' endstr]); - if abs(zshift)>1e-5 - % calculate new psiscatvol - psitdi=tdi('\results::psi'); - rmesh=psitdi.dim{1}; - zmesh=psitdi.dim{2}; - zthom=mdsdata('dim_of(\thomson:te,1)'); - zeffshift=zshift; - % set zeffshift time array same as psitdi - switch length(zeffshift) - case 1 - zeffshift=zeffshift * ones(size(psitdi.dim{3})); - case length(psitdi.dim{3}) - % ok - case length(psiscatvol.dim{1}) - zeffshift=interp1(psiscatvol.dim{1},zeffshift,psitdi.dim{3}); - otherwise - disp(' bad time dimension for zshift') - disp(['it should be 1 or ' num2str(length(psiscatvol.dim{1})) ' or ' num2str(length(psitdi.dim{3}))]) - end - for it=1:length(psiscatvol.dim{1}) - itpsitdi=iround(psitdi.dim{3},psiscatvol.dim{1}(it)); - psirz=psitdi.data(:,:,itpsitdi); - psiscatvol0=griddata(rmesh,zmesh,psirz',0.9*ones(size(zthom)),zthom-zeffshift(itpsitdi)); - psiscatvol.data(it,:)=psiscatvol0; - end - end - if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data) - for ir=1:length(psiscatvol.dim{2}) - rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))'; - end - else - rho=NaN; - end - trace.dim=[{rho};{time}]; - trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}]; - trace.x=rho; - trace.t=time; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - mdsclose; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'profnerho','profterho'} - % vol from psitbx - mdsopen(shot); - error_status=1; - if strcmp(TCVkeywrdcase{index},'profnerho') - nodenameeff=['\results::THOMSON.PROFILES.AUTO:ne']; - avers=tdi('\results::THOMSON.PROFILES.AUTO:ne:version_num'); - end - if strcmp(TCVkeywrdcase{index},'profterho') - nodenameeff=['\results::THOMSON.PROFILES.AUTO:te']; - avers=tdi('\results::THOMSON.PROFILES.AUTO:te:version_num'); - end - if avers.data==0 - % may be because nodes not yet filled in, so call once a node - ab=tdi(nodenameeff); - avers=tdi([nodenameeff ':version_num']); - end - if avers.data>0 - tracetdi=tdi(nodenameeff); - if avers.data < 2.99 - % for earlier version the bug made it to have logically (rho,t) - if ~isempty(tracetdi.dim) && ~ischar(tracetdi.data) - trace.x=tracetdi.dim{1}; - trace.t=tracetdi.dim{2}; - error_status=0; - else - error_status=2; - trace.x=[]; - trace.t=[]; - end - else - trace.data=tracetdi.data'; % error in dimensions for autofits - if ~isempty(tracetdi.dim) && ~ischar(tracetdi.data) - disp('assumes dim{2} for x in THOMSON.PROFILES.AUTO') - trace.x=tracetdi.dim{2}; - trace.t=tracetdi.dim{1}; - error_status=0; - else - trace.x=[]; - trace.t=[]; - error_status=2; - end - end - else - tracetdi=avers; - trace.x=[]; - trace.t=[]; - end - trace.dim=[{trace.x};{trace.t}]; - trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}]; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) '; Thomson autfits from ' nodenameeff ]; - mdsclose; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'B0'} - % B0 at R0=0.88 - if liuqe_version==-1 - % mdsopen( 'pcs', shot); %synthetic shot generated with FBTE and MGAMS. - mdsopen(shot) - nodenameeff = 'tcv_eq("BZERO","FBTE")'; - tracetdi=tdi(nodenameeff); - trace.data = tracetdi.data; - trace.t = tracetdi.dim{1}; - - else - mdsopen(shot); - nodenameeff=['\magnetics::iphi']; - tracetdi=tdi(nodenameeff); - R0EXP=0.88; - trace.data=192.E-07 * 0.996 *tracetdi.data/R0EXP; - trace.t=tracetdi.dim{1}; - end - trace.x=[]; -%keyboard - trace.dim=[{trace.t}]; - trace.dimunits=[{'time [s]'}]; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) ';' nodenameeff]; - mdsclose; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'qrho'} - % q profile on psi from liuqe - if liuqe_version==-1 - mdsopen( 'pcs', shot); %synthetic shot generated with FBTE and MGAMS. - else - mdsopen(shot); - end - nodenameeff=[begstr 'q_psi' endstr]; - tracetdi=tdi(nodenameeff); - trace.data=tracetdi.data; - trace.x=sqrt(tracetdi.dim{1}/(length(tracetdi.dim{1})-1)); - trace.t=tracetdi.dim{2}; - trace.dim=[{trace.x};{trace.t}]; - trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}]; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) ';' nodenameeff]; - mdsclose; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'vol'} - % vol from psitbx - if liuqe_version==-1 - begstr = 'tcv_psitbx("'; - mdsopen( 'pcs', shot); %synthetic shot generated with FBTE and MGAMS. - nodenameeff=[begstr 'vol' endstr]; - tracetdi=tdi(nodenameeff); - else - mdsopen(shot); - nodenameeff=['\results::psitbx:vol']; - tracetdi=tdi(nodenameeff); - if i_23==2 - ['LIUQE' endstr(2:2)] - psi=psitbxtcv(shot,'+0',['LIUQE' endstr(2:2)]); - fsd=psitbxp2p(psi,'FS'); - fsg=psitbxfsg(fsd); - tracetdi.data = fsg.vol.x; - tracetdi.dim{1}=fsg.vol.grid.x{:}'; - tracetdi.dim{2}=fsg.vol.grid.t'; - end - end - trace.data=tracetdi.data; - if isempty(tracetdi.data) || ischar(tracetdi.data) - trace.x=tracetdi.dim; - trace.t=tracetdi.dim; - trace.dim=tracetdi.dim; - else - trace.x=tracetdi.dim{1}; - trace.t=tracetdi.dim{2}; - trace.dim=[{trace.x};{trace.t}]; - end - trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}]; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) ';' nodenameeff '_' num2str(liuqe_version)]; - mdsclose; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'rhovol'} - % vol from psitbx - if liuqe_version==-1 - begstr = 'tcv_psitbx("'; - mdsopen( 'pcs', shot); %synthetic shot generated with FBTE and MGAMS. - nodenameeff=[begstr 'vol' endstr]; - tracetdi=tdi(nodenameeff); - else - mdsopen(shot); - nodenameeff=['\results::psitbx:vol']; - tracetdi=tdi(nodenameeff); - if i_23==2 - ['LIUQE' endstr(2:2)] - psi=psitbxtcv(shot,'+0',['LIUQE' endstr(2:2)]); - fsd=psitbxp2p(psi,'FS'); - fsg=psitbxfsg(fsd); - tracetdi.data = fsg.vol.x; - tracetdi.dim{1}=fsg.vol.grid.x{:}'; - tracetdi.dim{2}=fsg.vol.grid.t'; - end - end - trace.data=tracetdi.data; - for i=1:size(tracetdi.data,2) - trace.data(:,i)=sqrt(tracetdi.data(:,i)./tracetdi.data(end,i)); - end - trace.x=tracetdi.dim{1}; - trace.t=tracetdi.dim{2}; - trace.dim=[{trace.x};{trace.t}]; - trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}]; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) '; sqrt(V/Va) from ' nodenameeff '_' num2str(liuqe_version)]; - mdsclose; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'sxr','sxR'} - % load TCV soft x-ray data - if nargin>=3 & ~isempty(varargin{1}) - i1=varargin{1}(1); - i2=varargin{1}(2); - else - i1=1; - i2=20; - end - if nargin>=4 & ~isempty(varargin{2}) - status=varargin{2}; - else - status=ones(i2-i1+1,1); - end - - % camera selection: 1-10. each camera has 20 channels - icamera=[0 1 0 0 0 0 0 0 0 0]; %index of the camera to use - if ~isempty(find(status(1:20*icamera*ones(10,1)) == 1)) - [fans,vangle,xchord,ychord,aomega,angfact]=xtomo_geometry(1,icamera); - % calculating intersection of the view lines with magnetic axis - if strcmp(data_type_eff,'sxR') - if nargin>=5 & ~isempty(varargin{3}) - zmag=varargin{3}; - else - zmag=loadTCVdata(shot,['zmag' '_' num2str(liuqe_version)]); - end - t_1=zmag.t(1); - t_2=zmag.t(end); - [xtomo_signal,t]=get_xtomo_data(shot,t_1,t_2,13e-6*16,icamera,angfact); - data=interp1(zmag.t,zmag.data,t'); - radius.data=VsxrTCVradius(data,xchord,ychord)'; - radius.t=t'; - varargout{1}={radius}; - trace.R=radius.data; - else - t_1=0.001; - t_2=3; - [xtomo_signal,t]=get_xtomo_data(shot,t_1,t_2,13e-6*16, ... - icamera,angfact); - end - end - trace.data=xtomo_signal; - trace.x=[1:size(trace.data,1)]'; - trace.t=t'; - trace.dim=[{trace.x} ; {trace.t}]; - trace.dimunits=[{'channel #'} ; {'time [s]'}]; - trace.name=[num2str(shot) ';' 'get_xtomo_data']; - error_status=0; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case 'ece' - % load TCV ECE data - % Status=1 => Not Read Yet - mdsopen(shot); - if ~isempty(find(status == 1)) - if eval(['~mdsdata(''node_exists("\\RESULTS::ECE:rho")'')']) - disp(['node \RESULTS::ECE:rho does not exist for shot = ' num2str(shot)]) - return - end - if eval(['mdsdata(''getnci("\\RESULTS::ECE:rho","length")'')==0']) - disp(['no data for \RESULTS::ECE:rho for shot = ' num2str(shot)]) - return - end - [TE_ECE,TE_ECE_ERR,RHO,R,T,TE_THOM,TE_THOM_ERR,Fcentral,CAL]=ece_te ... - (shot,[0.1 0.29],10,10); - end - a=min(find(R(:,1)>=0)); - b=max(find(R(:,1)>=0)); - trace.data=TE_ECE(a:b,:)'; - trace.t=T(a:b); - trace.x=[1:size(trace.data,1)]'; - trace.dim=[{trace.x} ; {trace.t}]; - trace.dimunits=[{'channel #'} ; {'time [s]'}]; - trace.R=R(a:b,:)'; - trace.name=[num2str(shot) ';' '\results::ece...']; - radius.data=trace.R; - radius.t=trace.t; - varargout{1}={radius}; - error_status=0; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - mdsclose; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case 'MPX' - % load TCV MPX data - % Status=1 => Not Read Yet - if isempty(zmag) - zmag=loadTCVdata(shot,'zmag'); - end - t_1=zmag.t(1); - t_2=zmag.t(end); - if ~isempty(find(status == 1)) - mdsopen(shot); - keyboard - signal=get_mds_mio('MPX',[t_1 t_2]); - mdsclose; - trace.data=signal.data; - for i=1:size(signal.dim{2},2) - trace.t(:,i)=signal.dim{1}; - end - end - trace.dim{1}={trace.t}; - trace.dimunits={'time [s]'}; - trace.name=[num2str(shot) ';' 'get_mds_mio(MPX)']; - [xchord,ychord]=mpx_geometry; - varargout{1}={VsxrTCVradius(zmag.data,xchord,ychord)}; - error_status=0; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case 'IOH' - % Ohmic transformer current - mdsopen(shot); - nodenameeff=[{'\magnetics::ipol[*,$1]'} {'OH_001'}]; - tracetdi=tdi(nodenameeff{:}); - trace.data=tracetdi.data; - trace.x=[]; - trace.t=tracetdi.dim{1}; - trace.dim=tracetdi.dim; - trace.dimunits={'time [s]'}; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) ';' nodenameeff{1} ',' nodenameeff{2}]; - mdsclose; - error_status=0; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case 'vloop' - % Loop voltage - mdsopen(shot); - nodenameeff=[{'\magnetics::vloop[*,$1]'} {'001'}]; - tracetdi=tdi(nodenameeff{:}); - trace.data=tracetdi.data; - trace.x=[]; - trace.t=tracetdi.dim{1}; - trace.dim=tracetdi.dim; - trace.dimunits={'time [s]'}; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) ';' nodenameeff{1} ',' nodenameeff{2}]; - mdsclose; - error_status=0; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case 'pgyro' - % ECH power for each gyro(1:9) and total (10) - mdsopen(shot); - nodenameeff=[{'\results::toray.input:p_gyro'}]; - tracetdi=tdi(nodenameeff{:}); - trace.data=tracetdi.data; - trace.x=[]; - trace.t=tracetdi.dim{1}; - trace.dim=tracetdi.dim; - trace.dimunits=tracetdi.dimunits; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) ';' nodenameeff]; - mdsclose; - error_status=0; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case 'jtor' - % \results::j_tor , 3-D so need to specify time dim index - mdsopen(shot); - nodenameeff=[{'\results::j_tor'} endstr]; - tracetdi=tdi(nodenameeff{:}); - trace.data=tracetdi.data; - trace.x=tracetdi.dim{1}; - trace.t=tracetdi.dim{3}; - trace.dim=tracetdi.dim; - trace.dimunits=tracetdi.dimunits; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) ';' nodenameeff]; - mdsclose; - error_status=0; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'neft:trial','teft:trial','neftav:trial','teftav:trial'} - % trial indx - mdsopen(shot); - eval(['nodenameeff={''' TCVsiglocation{index} ':trial''};']); - tracetdi=tdi(nodenameeff{:}); - if isempty(trialindx) - error('trialindx should not be empty, check call or ask O. Sauter'); - end - trace.data=tracetdi.data(:,:,trialindx+1); - trace.x=tracetdi.dim{1}; - trace.t=tracetdi.dim{2}; - trace.dim=tracetdi.dim(1:2); - trace.dimunits=tracetdi.dimunits(1:2); - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) ' ; ' nodenameeff{:} ' ; trialindx=' num2str(trialindx) ]; - mdsclose; - error_status=0; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - case {'vi_tor', 'vi_torfit', 'vi_pol', 'vi_polfit', 'Ti', 'Tifit', 'ni', 'nifit', 'zeffcxrs', 'zeffcxrsfit'} - proffit = ''; - kwd_eff = TCVkeywrdcase{index}; - ii=strfind(kwd_eff,'fit'); - if ~isempty(ii); - proffit = ':proffit'; - kwd_eff = kwd_eff(1:ii-1); - end - if strcmp(kwd_eff,'zeffcxrs'); kwd_eff = 'zeff'; end - mdsopen(shot); - eval(['nodenameeff=''\results::cxrs' proffit ':' kwd_eff endstr ''';']); - tracetdi=tdi(nodenameeff); - trace.data=tracetdi.data; - if length(tracetdi.dim)>1; - trace.x=tracetdi.dim{1}; - trace.t=tracetdi.dim{2}; - else - trace.x=[]; - trace.t=[]; - end - trace.dim=tracetdi.dim; - trace.dimunits=tracetdi.dimunits; - % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus isfield using isa does not work - if any(strcmp(fieldnames(tracetdi),'units')) - trace.units=tracetdi.units; - end - trace.name=[num2str(shot) ';' nodenameeff]; - % add error bars - eval(['nodenameeff=''\results::cxrs' proffit ':' kwd_eff ':err' endstr ''';']); - nodenameeff=[{'\results::cxrs:vi_tor:err'} endstr]; - tracetdi=tdi(nodenameeff{:}); - trace.std = tracetdi.data; - if length(tracetdi.dim)>2; - trace.std_rho = tracetdi.dim{1}; - trace.std_t = tracetdi.dim{2}; - else - trace.std_rho = []; - trace.std_t = []; - end - mdsclose; - error_status=0; - - %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& - otherwise - % eval(['!mailto_Andrea ''from loadTCVdata, data_type_eff= ' data_type_eff '''']) - disp(['this data_type_eff' ' ' data_type_eff ' ' 'not yet programmed in loadTCVdata, ask Andrea.Scarabosio@epfl.ch']); - -end diff --git a/matlab/TCV/tcv_requests_mapping.m b/matlab/TCV/tcv_requests_mapping.m index 9909a953e12faf245943034c299fe7db3eee168c..47be410cad59e6f6fdceca56a4e9ae9c05e55821 100644 --- a/matlab/TCV/tcv_requests_mapping.m +++ b/matlab/TCV/tcv_requests_mapping.m @@ -457,6 +457,11 @@ switch lower(data_request) 'bb=interpos(aa.t,aa.data(ix,:),gdat_tmp.t,-100);' ... 'gdat_tmp.data = 4e-7*pi.*gdat_tmp.data.^2.*bb./1.22; gdat_tmp.label=''\tau_R'';gdat_tmp.gdat_request=''tau_r'';' ... 'gdat_tmp.units=''s'';gdat_tmp.help=''tauR=mu0 a^2 / 1.22 etaneo(rhopol=0.4)''']; + case 'tau_tot' + mapping.timedim = 1; + mapping.label = 'total energy confinment time'; + mapping.method = 'tdi'; + mapping.expression = '\results::conf:tau'; case 'te' mapping.timedim = 2; mapping.label = 'Te'; diff --git a/matlab/TCV/xtomo_geometry.m b/matlab/TCV/xtomo_geometry.m deleted file mode 100644 index a809b9e65cca29d8df8eff407eaac96d7b7ebc4c..0000000000000000000000000000000000000000 --- a/matlab/TCV/xtomo_geometry.m +++ /dev/null @@ -1,402 +0,0 @@ -function [fans,vangle,xchord,ychord,aomega,angfact]=xtomo_geometry(i_detec,fans) - -% ----[anton.public] -% -% function -% -% [fans,vangle,xchord,ychord,aomega,angfact]=xtomo_geometry(i_detec,fans); -% -% inputs: -% -% i_detec: =2: Xtomo prototype cameras (shot# < 6768) -% =1: Xtomo 9-cameras (shot# > 682x) -% -% outputs: -% -% fans: camera switch, 1=on,0=off (1x10) -% vangle: angle between detect. surface normal and pos. x-axis (1x10) -% xchord: two x-coordinates (2xnl) and -% ychord: two y-coord. for each line (2xnl), they specify start + end points -% aomega: etendue in mm^2 x steradians -% angfact: angular factors, inverse of relative etendue (throughput) (20x10) -% -% uses: -% AOMEGA=etendue_n2(b1x,b1y,b1z,b2x,b2y,b2z,z01,z02,X0,cw); -% angular_fact_*.mat , '*'=i_detec -% -% -%---------------- M.Anton 14/3/95 ------------------------------------------- - -disp('*----------------------------*') -disp('| this is xtomo_geometry |') -disp('*----------------------------*') - -global xap yap xdet ydet -global ae da -% ======== tokamak parameters ================================================ - -load tcv_vesc1 - - - - -xcont=rzvin(:,1); -ycont=rzvin(:,2); -xmin=min(xcont); -ymin=min(ycont); -xmax=max(xcont); -ymax=max(ycont); -xedge=100; -yedge=60; - - - - -% ========= detector parameters ============================================= - -if i_detec==2 - - cw=1; % detector numbers cw=1:clockwise cw=0:ccw - - if nargin<2 - fans=[0 0 0 0 0 0 1 0 1 0]; % camera switch - end - - vangle=[90 90 90 0 0 0 0 -90 -90 -90]; - % angle of detector surface normal - xpos=[0 0 0 0 0 0 118.05 0 87.84 0]; - % x position of the diaphragmas in [cm] - ypos=[0 0 0 0 0 0 -46 0 -80.45 0]; - % y position of the diaphragmas in [cm] - ae=[0 0 0 0 0 0 -2.5 0 -0.1 0]/10; - %excentricity array/diaphragm in [cm] - da=[0 0 0 0 0 0 10.1 0 -11.7 0]/10; - % diaphragma-array distance in [cm] - da2=da; - - d1=0.950; % detector width in mm - d2=4.050; % detector length in mm - b1=1.000; % aperture width in mm - b2=4.000*ones(1,10); % aperture length in mm - b3x=0; % aperture thickness in mm - b3y=0; - -elseif i_detec==1 - - cw=1; % detector numbers cw=1:clockwise cw=0:ccw - - if nargin<2 - fans=[1 1 1 1 1 1 1 1 0 1]; % camera switch - end - - vangle=[90 90 90 0 0 0 0 -90 -90 -90]; - % angle of detector surface normal - xpos=[71.5 87.7 103.9 123.1 123.1 123.1 123.1 104.04 87.84 71.64]; - % x position of the diaphragmas in [cm] - ypos=[80.35 80.35 80.35 48.1 1.95 -2.45 -48.6 -80.35 -80.35 -80.35]; - % y position of the diaphragmas in [cm] - ae=[-8 0 8 5 9 -9 -5 8 0 -8]/10; %excentricity array/diaphragm [cm] - - - ae= ae + [-0.0915 0 0.1361 0.2123 0.0923 -0.0994 ... - 0.0872 -0.1520 0 0.9410 ]/10; - ae(1)=ae(1)+0.1/10; - ae(3)=8/10+0.14/10; - ae(4)=4.9/10; - ae(5)=9/10+0.2/10; - ae(6)=ae(6)-0.2/10; - ae(7)=-4.9/10; - ae(10)=-7.1/10; - - da= [12.4 9.9 12.4 9.9 13.4 13.4 9.9 -12.4 -9.9 -12.4]/10; - % diaphragma-array distance in [cm] (poloidal) - da2=[37 34.4 37 55.9 59.4 59.4 55.9 37 34.4 37]/10; - % dist to diaphragm in toroidal direction [cm]; - deltada=[ -0.0311 0 -0.0458 -0.1179 -0.0615 -0.1105 ... - -0.0510 -0.0515 0 -0.3223]/10; - deltada(4)=0; - deltada(6)=0; - - - da=da+ deltada; - - da2=da2+deltada; - - - d1=0.90; % detector width in mm - d2=4.0; % detector length in mm - b1=0.800; % aperture width in mm (pol.) - b2=[8 8 8 15 15 15 15 8 8 8]; - % aperture length in mm (tor.) - b3x=0.020; % aperture thickness in mm (poloidal) - b3y=0; % aperture thickness in mm (toroidal) -end - - - -%======== calculation of the chords of view =================================== - -nact=sum(fans); -iact=find(fans); -ndet=20; -ncam=10; - - - -% ---- apertures: ------------------ - -xap=ones(ndet,1)*xpos(iact); -xap=xap(:)'; -yap=ones(ndet,1)*ypos(iact); -yap=yap(:)'; - -% ---- detectors: ------------------ - -vorz(find(vangle==90))=(-1)^(cw+1)*ones(size(find(vangle==90))); -vorz(find(vangle==0))=(-1)^cw*ones(size(find(vangle==0))); -vorz(find(vangle==-90))=(-1)^cw*ones(size(find(vangle==-90))); - - -dete=(-9.025:0.950:9.025)'/10*vorz(iact)+ones(ndet,1)*ae(iact); - -dum_ae=dete(:)'; - -dum_vangle=ones(ndet,1)*vangle(iact); -dum_vangle=dum_vangle(:)'; - - - -ivert=find(dum_vangle==90 | dum_vangle==-90); -ihori=find(dum_vangle==0); - -dum_da=ones(ndet,1)*da(iact); -dum_da=dum_da(:)'; - -dxd=zeros(1,ndet*nact); -dyd=zeros(1,ndet*nact); - -dxd(ivert)=dum_ae(ivert); -dxd(ihori)=dum_da(ihori); - -dyd(ivert)=dum_da(ivert); -dyd(ihori)=dum_ae(ihori); - -xdet=xap+dxd; -ydet=yap+dyd; - - -%plot_vessel(rzvin,rzvout) -%hold on -% plot(xap,yap,'.g',xdet,ydet,'.m') - - -% ---- calculate the equations of lines of sight - -m=(ydet-yap)./(xdet-xap); -b=ydet-m.*xdet; - -nl=length(xdet); -xchord=zeros(2,nl); -ychord=zeros(2,nl); - - -xchord(1,:)=xdet;ychord(1,:)=ydet; - -iup=find(dum_vangle==90); -isi=find(dum_vangle==0); -ido=find(dum_vangle==-90); - -if ~isempty(iup) -ychord(2,iup)=ymin*ones(size(iup)); -xchord(2,iup)=(ychord(2,iup)-b(iup))./m(iup); -end -if ~isempty(ido) -ychord(2,ido)=ymax*ones(size(ido)); -xchord(2,ido)=(ychord(2,ido)-b(ido))./m(ido); -end -if ~isempty(isi) -xchord(2,isi)=xmin*ones(size(isi)); -ychord(2,isi)=m(isi).*xchord(2,isi)+b(isi); -end - -ileft=find(xchord(2,:)<xmin); - -if ~isempty(ileft) -xchord(2,ileft)=xmin*ones(size(ileft)); -ychord(2,ileft)=m(ileft).*xchord(2,ileft)+b(ileft); -end -irig=find(xchord(2,:)>xmax); -if ~isempty(irig) -xchord(2,irig)=xmax*ones(size(irig)); -ychord(2,irig)=m(irig).*xchord(2,irig)+b(irig); -end - - - -%======== prepare output ====================================================== - -vangle=vangle(iact); - -%======== calculation of angular correction factors, if necessary ============= - -if i_detec==2 & exist('angular_fact_2.mat')==2 - - disp('loading angular_fact_2') - load angular_fact_2 - - -elseif i_detec==1 & exist('angular_fact_1.mat')==2 - - disp('loading angular_fact_1') - load angular_fact_1 - -else - - aomega=zeros(ndet,ncam); - angfact=ones(ndet,ncam); - - - for l=1:sum(fans) - -% Z0X=abs(da(iact(l))*10) -% Z0Y=abs(da2(iact(l))*10) -% X0=ae(iact(l))*10 % back to mm, sorry about that... -% X0=X0*vorz(iact(l)) -% B2=b2(iact(l)) -% AOMEGA=etendue_n(b1,B2,b3x,b3y,Z0X,Z0Y,X0,cw); - - b1x=0.8; - b1y=6; - b1z=0.02; - b2x=10000000; - b2y=b2(iact(l)); - b2z=0; - z01=abs(da(iact(l))*10); - z02=abs(da2(iact(l))*10); - X0=ae(iact(l))*10*vorz(iact(l)); -keyboard - AOMEGA=etendue_n2(b1x,b1y,b1z,b2x,b2y,b2z,z01,z02,X0,cw); - - aomega(:,iact(l))=AOMEGA(:,1); - - - end - - - indm=min(find(aomega==max(aomega(:)))); - aomegan=aomega/aomega(indm); - nonz=find(aomega); - angfact(nonz)=ones(size(nonz))./aomegan(nonz); - - angfact=round(1000*angfact)/1000; - - - unitstring='units aomega: mm^2 * sterad'; - - if i_detec==1 - save angular_fact_1 angfact aomega unitstring - elseif i_detec==2 - save angular_fact_2 angfact aomega unitstring - end - -end -return - -th=atan(diff(ychord)./diff(xchord)); -thp=th; -neg=find(thp<0); -thp(neg)=180+thp(neg); -thdet=ones(1,20*sum(fans)); -for k=1:sum(fans) - thdet((k-1)*20+1:k*20)=vangle(k)*ones(1,20); -end -angles=thdet-thp; -mist=find(angles<0 & abs(angles)>90); -angles(mist)=angles(mist)+180; -th_inc=angles*pi/180; - - -% ---- correct for the edges of tcv ( some chords may be too long ) - - - - - - down=find(xcont>xedge & ycont<-yedge); - up=find(xcont>xedge & ycont>yedge); - cd=polyfit(xcont(down),ycont(down),1); - cu=polyfit(xcont(up),ycont(up),1); - - - iu1=find(xchord(1,:)>xedge & ychord(1,:)>0 & dum_vangle==-90 ); - if ~isempty(iu1) - xchord(1,iu1)=-(b(iu1)-cu(2))./(m(iu1)-cu(1)+eps); - ychord(1,iu1)=m(iu1).*xchord(1,iu1)+b(iu1); - end - - iu2=find(xchord(2,:)>xedge & ychord(2,:)>0 & ychord(1,:) & .... - dum_vangle==-90); - if ~isempty(iu2) - xchord(2,iu2)=-(b(iu2)-cu(2))./(m(iu2)-cu(1)+eps); - ychord(2,iu2)=m(iu2).*xchord(2,iu2)+b(iu2); - end - - id1=find(xchord(1,:)>xedge & ychord(1,:)<0 & dum_vangle==90); - if ~isempty(id1) - xchord(1,id1)=-(b(id1)-cd(2))./(m(id1)-cd(1)+eps); - ychord(1,id1)=m(id1).*xchord(1,id1)+b(id1); - end - - id2=find(xchord(2,:)>xedge & ychord(2,:)<0 & dum_vangle==90); - if ~isempty(id2) - xchord(2,id2)=-(b(id2)-cd(2))./(m(id2)-cd(1)+eps); - ychord(2,id2)=m(id2).*xchord(2,id2)+b(id2); - end - - ilow=find(ychord(1,:)<ymin); - ihig=find(ychord(1,:)>ymax); - ilef=find(xchord(1,:)<xmin); - irig=find(xchord(1,:)>xmax); - if ~isempty(ilow) - ychord(1,ilow)=ymin*ones(size(ilow)); - xchord(1,ilow)=ymin./m(ilow)-b(ilow)./m(ilow); - end - if ~isempty(ihig) - ychord(1,ihig)=ymax*ones(size(ihig)); - xchord(1,ihig)=ymax./m(ihig)-b(ihig)./m(ihig); - end - if ~isempty(ilef) - xchord(1,ilef)=xmin*ones(size(ilef)); - ychord(1,ilef)=m(ilef)*xmin+b(ilef); - end - if ~isempty(irig) - xchord(1,irig)=xmax*ones(size(irig)); - ychord(1,irig)=m(irig)*xmax+b(irig); - end - - - ilow=find(ychord(2,:)<ymin); - ihig=find(ychord(2,:)>ymax); - ilef=find(xchord(2,:)<xmin); - irig=find(xchord(2,:)>xmax); - if ~isempty(ilow) - ychord(2,ilow)=ymin*ones(size(ilow)); - xchord(2,ilow)=ymin./m(ilow)-b(ilow)./m(ilow); - end - if ~isempty(ihig) - ychord(2,ihig)=ymax*ones(size(ihig)); - xchord(2,ihig)=ymax./m(ihig)-b(ihig)./m(ihig); - end - if ~isempty(ilef) - xchord(2,ilef)=xmin*ones(size(ilef)); - ychord(2,ilef)=m(ilef)*xmin+b(ilef); - end - if ~isempty(irig) - xchord(2,irig)=xmax*ones(size(irig)); - ychord(2,irig)=m(irig)*xmax+b(irig); - end - - - - - diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m new file mode 100644 index 0000000000000000000000000000000000000000..710b3e90a941e3f9104a60862a60ddf0f3a8a4f6 --- /dev/null +++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m @@ -0,0 +1,149 @@ +function [ids_core_sources,ids_core_sources_description,varargout] = ... + tcv_get_ids_core_sources(shot,ids_equil_empty,gdat_params,varargin) +% +% [ids_core_sources,ids_core_sources_description,varargout] = ... +% tcv_get_ids_core_sources(shot,ids_equil_empty,gdat_params,varargin); +% +% +% gdat_params: gdat_data.gdat_params to get all params passed from original call +% + +machine = 'tcv'; +tens_time = -1; +tens_rho = -0.1; + +if exist('gdat_params','var') + [ids_core_sources, params_core_sources] = ... + tcv_ids_headpart(shot,ids_equil_empty,'cores_sources','gdat_params',gdat_params,varargin{:}); +else + [ids_core_sources, params_core_sources] = ... + tcv_ids_headpart(shot,ids_equil_empty,'cores_sources',varargin{:}); + aa=gdat_tcv; + gdat_params = aa.gdat_params; % to get default params +end +params_eff_ref = gdat_params; params_eff_ref.doplot=0; +try params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def) + +% initialize description +ids_core_sources_description = []; + +%% +% save source template to fill for differnet source options +source_template = ids_core_sources.source{1}; +profiles_template = source_template.profiles_1d{1}; +globals_template = source_template.global_quantities{1}; +%% get time grid and data form gdat & initialize source from template +ohm_gdat = gdat(shot,'ohm_data'); +% fix time grid to ohm time grid, same as bs +t_grid = ohm_gdat.t; n_t = numel(t_grid); t_spacing = mean(diff(t_grid)); + +% ohm +id_ohm.description = 'Source from ohmic heating'; +id_ohm.index = 7; id_ohm.name = 'ohmic'; +ids_core_sources.source{7} = source_template; +ids_core_sources.source{7}.identifier = id_ohm; +ids_core_sources.source{7}.profiles_1d(1:n_t) = {profiles_template}; +ids_core_sources.source{7}.global_quantities(1:n_t) = {globals_template}; +ohm_data = ohm_gdat.ohm.ohm_data; + +% ec +id_ec.description = 'Sources from electron cyclotron heating and current drive'; +id_ec.index = 3; id_ec.name = 'ec'; +ids_core_sources.source{3} = source_template; +ids_core_sources.source{3}.identifier = id_ec; +ids_core_sources.source{3}.profiles_1d(1:n_t) = {profiles_template}; +ids_core_sources.source{3}.global_quantities(1:n_t) = {globals_template}; +% load data +ec_gdat = gdat(shot,'ec_data'); ec_data = ec_gdat.ec.ec_data; +ec_time = ec_gdat.ec.ec_data.cd_dens.t; +% interpolate totals on t_grid +ec_total_pow = interp1(ec_data.p_abs_plasma.t,ec_data.p_abs_plasma.data(12,:),t_grid); +ec_total_pow(isnan(ec_total_pow)) = 0; +ec_total_cur = interp1(ec_data.p_abs_plasma.t,ec_data.cd_tot.data(12,:),t_grid); +ec_total_cur(isnan(ec_total_cur)) = 0; + +% bs +id_bs.description = 'Bootstrap current'; +id_bs.index = 13; id_bs.name = 'bootstrap_cuurent'; +ids_core_sources.source{13} = source_template; +ids_core_sources.source{13}.identifier = id_bs; +ids_core_sources.source{13}.profiles_1d(1:n_t) = {profiles_template}; +ids_core_sources.source{13}.global_quantities(1:n_t) = {globals_template}; +% load data +bs_gdat = gdat(shot,'bs_data'); bs_data = bs_gdat.bs.bs_data; + +% nbi +ids_core_sources.source{2} = source_template; +% to be filled + +% total +ids_core_sources.source{1} = source_template; +% to be filled +powers_gdat = gdat(shot,'powers'); + +%% fill profiles for times n t_grid +for ii = 1:n_t + % ohm + % 1d_profiles + ids_core_sources.source{7}.profiles_1d{ii}.time = t_grid(ii); + ids_core_sources.source{7}.profiles_1d{ii}.j_parallel = ... + ohm_data.cd_dens.data(:,ii)'; + % globals + ids_core_sources.source{7}.global_quantities{ii}.time = t_grid(ii); + ids_core_sources.source{7}.global_quantities{ii}.power = powers_gdat.ohm.data(ii); + + % ec, need to find matching time index\ + ids_core_sources.source{3}.profiles_1d{ii}.time = t_grid(ii); % profiles time + ids_core_sources.source{3}.global_quantities{ii}.time = t_grid(ii); % globals time + [~, ec_idx] = min(abs(ec_time-t_grid(ii))); + if ec_time(ec_idx) == t_grid(ii) % same grid + % current density + ids_core_sources.source{3}.profiles_1d{ii}.j_parallel = ... + ec_data.cd_dens.data(:,end,ec_idx)'; + % integrated current density + ids_core_sources.source{3}.profiles_1d{ii}.current_parallel_inside = ... + ec_data.cd_integrated.data(:,ii)'; + + else + ids_core_sources.source{3}.profiles_1d{ii}.j_parallel = ... + zeros(1,numel(ec_data.cd_dens.x(1,:)))'; + ids_core_sources.source{3}.profiles_1d{ii}.current_parallel_inside = ... + zeros(1,numel(ec_data.cd_integrated.x(1,:)))'; + end + % globals + ids_core_sources.source{3}.global_quantities{ii}.time = t_grid(ii); + ids_core_sources.source{3}.global_quantities{ii}.power = ec_total_pow(ii); + ids_core_sources.source{3}.global_quantities{ii}.current_parallel = ec_total_cur(ii); + + % bs + ids_core_sources.source{13}.profiles_1d{ii}.j_parallel = ... + bs_data.cd_dens.data(:,ii)'; + +end + +%% add descriptions for profiles_1d + +desc.source = 'available by now, 1:total, 3:ec, 7:ohm, 13:bootstrap'; +desc.profiles_1d.j_parallel = 'parallel current density'; +desc.profiles_1d.current_parallel_inside = 'integrated parallel current density'; +desc.globals.power = 'total power coupled to the plasma'; +desc.globals.current_parallel = 'total parallel current driven'; +ids_core_sources_description = desc; + +%% +if nargin <= 2 + ids_core_sources.code.name = ['tcv_get_ids_core_sources, within gdat, with shot= ' num2str(params_core_sources.shot) ]; +else + ids_core_sources.code.name = ['tcv_get_ids_core_sources, within gdat, with shot= ' num2str(params_core_sources.shot) '; varargin: ' varargin{:}]; +end +ids_core_sources_description.code.name = ids_core_sources.code.name; + +ids_core_sources.code.output_flag = zeros(size(ids_core_sources.time)); + +% cocos automatic transform +if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic')) + [ids_core_sources,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_core_sources,'core_sources',gdat_params.cocos_in, ... + gdat_params.cocos_out,gdat_params.ipsign_out,gdat_params.b0sign_out,gdat_params.ipsign_in,gdat_params.b0sign_in, ... + gdat_params.error_bar,gdat_params.nverbose); + +end \ No newline at end of file diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_transport.m b/matlab/TCV_IMAS/tcv_get_ids_core_transport.m index 48edb64fd64939992aadf5ddcbefca18af3ebc33..8d0522c92c3fb713c926d294cbefddce64b7c283 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_core_transport.m +++ b/matlab/TCV_IMAS/tcv_get_ids_core_transport.m @@ -1,8 +1,8 @@ function [ids_core_transport,ids_core_transport_description,varargout] = ... tcv_get_ids_core_transport(shot,ids_equil_empty,gdat_params,varargin) % -% [ids_core_profiles,ids_core_profiles_description,varargout] = ... -% tcv_get_ids_core_profiles(shot,ids_equil_empty,gdat_params,varargin); +% [ids_core_transport,ids_core_transport_description,varargout] = ... +% tcv_get_ids_core_transport(shot,ids_equil_empty,gdat_params,varargin); % % % gdat_params: gdat_data.gdat_params to get all params passed from original call @@ -40,7 +40,7 @@ ids_core_transport.model{3}.name = comment; ids_core_transport.model{3}.identifier.index = 2; ids_core_transport.model{3}.name = 'transport_solver'; -% read data and setup time +% read data and setup time params_eff.data_request='\results::conf:chie'; temp_1d.chie = gdat(params_core_transport.shot,params_eff.data_request); temp_1d_desc.chie = params_eff.data_request; @@ -83,4 +83,4 @@ if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic')) gdat_params.cocos_out,gdat_params.ipsign_out,gdat_params.b0sign_out,gdat_params.ipsign_in,gdat_params.b0sign_in, ... gdat_params.error_bar,gdat_params.nverbose); -end \ No newline at end of file +end diff --git a/matlab/TCV_IMAS/tcv_get_ids_summary.m b/matlab/TCV_IMAS/tcv_get_ids_summary.m index 777eed4318e4fddfd68f5235a3880cf9e34b1d9b..06fb3c6a2dc9b17ad0ad50a94da1011673e07fff 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_summary.m +++ b/matlab/TCV_IMAS/tcv_get_ids_summary.m @@ -36,6 +36,7 @@ ids_summary.code.name = 'gdat'; ids_summary.code.version = git_release_hash; ids_summary.code.repository = git_url; +%% boundary % use equilibrium time as reference, so load first elongation(time) params_eff = params_eff_ref; @@ -59,7 +60,6 @@ ids_summary.time = unique([t_before; [kappa.t(1):max((kappa.t(end)-kappa.t(1))/( ids_summary_description.time = 'time from gdat(shot,''kappa'') and adding from 0 up to t(ip_trapeze(end))'; params_eff = params_eff_ref; -% boundary: boundary.elongation = kappa; boundary_desc.elongation = kappa.gdat_params.data_request; params_eff.data_request = 'r_geom'; @@ -121,13 +121,13 @@ if ~isempty(special_fields) % $$$ end end -% disruption, use specific function +%% disruption, use specific function % to be decided if changes only ids_summary_out.disruption or also elsewhere..., if ids_summary_out.time is changed, % then all the other data should be related, thus call this first? [ids_summary_out,ids_summary_out_description] = ids_summary_disruption(shot, ids_summary, gdat_params); -% gas: bottom/top?, deuterium, helium_xx, ? +%% gas: bottom/top?, deuterium, helium_xx, ? params_eff.data_request = 'gas_flux'; gas_injection_rates.deuterium = gdat(params_summary.shot,params_eff); gas_injection_rates_desc.deuterium = params_eff.data_request; @@ -148,7 +148,8 @@ for i=1:numel(gas_injection_rates_fieldnames) end end %params_eff.doplot=1; -% global_quantities: + +%% global_quantities params_eff.data_request = 'b0'; global_quantities.b0 = gdat(params_summary.shot,params_eff); % do not add 'source','liuqe' to get all times global_quantities_desc.b0 = [params_eff.data_request ' ; ' global_quantities.b0.data_fullpath]; @@ -182,6 +183,9 @@ global_quantities_desc.current_non_inductive = [params_eff.data_request ' ; ' gl params_eff.data_request = 'w_mhd'; global_quantities.energy_diamagnetic = gdat(params_summary.shot,params_eff); global_quantities_desc.energy_diamagnetic = [params_eff.data_request ' ; ' global_quantities.energy_diamagnetic.data_fullpath]; +params_eff.data_request = 'tau_tot'; +global_quantities.tau_energy = gdat(params_summary.shot,params_eff); +global_quantities_desc.tau_energy = [params_eff.data_request ' ; ' global_quantities.tau_energy.data_fullpath]; params_eff.data_request = '\tcv_shot::top.results.conf:we'; global_quantities.energy_electrons_thermal = gdat(params_summary.shot,params_eff); global_quantities_desc.energy_electrons_thermal = [params_eff.data_request ' ; ' global_quantities.energy_electrons_thermal.data_fullpath]; @@ -263,6 +267,7 @@ end ids_summary.global_quantities.r0.value = global_quantities.r0.data; ids_summary.global_quantities.r0.source = ['gdat request: ' global_quantities_desc.r0]; +%% volume average % cst zeff params_eff.data_request = '\tcv_shot::top.results.ibs:z_eff'; volume_average.zeff = gdat(params_summary.shot,params_eff); @@ -283,6 +288,29 @@ for i=1:numel(volume_average_fieldnames) end end +%% line average +params_eff.data_request = 'nel'; +line_average.n_e = gdat(params_summary.shot,params_eff); +line_average_desc.n_e = [params_eff.data_request ' ; ' volume_average.zeff.data_fullpath ' (not vol avrg but to match Iohm)']; + +special_fields = {}; % fields needing non-automatic treatments +line_average_fieldnames = fieldnames(line_average); +for i=1:numel(line_average_fieldnames) + if ~any(strcmp(line_average_fieldnames{i},special_fields)) + if ~isstruct(ids_summary.line_average.(line_average_fieldnames{i}).value) + ids_summary.line_average.(line_average_fieldnames{i}).value = interp1( ... + line_average.(line_average_fieldnames{i}).t,line_average.(line_average_fieldnames{i}).data, ... + ids_summary.time,'linear',NaN); + ids_summary.line_average.(line_average_fieldnames{i}).source = ['gdat request: ' line_average_desc.(line_average_fieldnames{i})]; + else + special_fields{end+1} = line_average_fieldnames{i}; + end + end +end + + +%% + % cocos automatic transform if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic')) [ids_summary,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_summary,'summary',gdat_params.cocos_in, ... diff --git a/matlab/TCV_IMAS/tcv_ids_headpart.m b/matlab/TCV_IMAS/tcv_ids_headpart.m index fa90b8cf836577146f499ba9cbe0ff7e2a53d311..53216596c61a9bbccaa9c08fc501e892ad01a1af 100644 --- a/matlab/TCV_IMAS/tcv_ids_headpart.m +++ b/matlab/TCV_IMAS/tcv_ids_headpart.m @@ -1,21 +1,20 @@ -function [ids_generic, params_ids_generic] = tcv_ids_headpart(shot,ids_in,ids_name,varargin); +function [ids_generic, params_ids_generic] = tcv_ids_headpart(shot,ids_in,ids_name,varargin) % % [ids_generic, params_ids_generic] = tcv_ids_headpart(shot,ids_in,ids_name,varargin); % % parses inputs and fill in ids_properties % -% % varargin options: -% -% 'comment': comment to include in ids_properties, using gdat_params for example cocos_in and cocos_out -% 'homogeneous_time': homogeneous_time in ids_properties: 1 (default) if the whole ids has same time, 0 otherwise -% 'gdat_params': gdat params structure +% 'comment': comment to include in ids_properties, using gdat_params +% for example cocos_in and cocos_out +% 'homogeneous_time': homogeneous_time in ids_properties: +% 1 (default) if the whole ids has same time, 0 otherwise +% 'gdat_params': gdat params structure % % % example: % [ids_equilibrium, params_ids_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium','comment','your comment'); % -% % initialize input parser p = inputParser;