diff --git a/crpptbx/TCV/loadTCVdata.m b/crpptbx/TCV/loadTCVdata.m index 16f90b2f0d59022174edb1f6bd70ba4cbcdb81be..43ed4a747f8dc20fe23ca46a3c9d4414740ad01c 100644 --- a/crpptbx/TCV/loadTCVdata.m +++ b/crpptbx/TCV/loadTCVdata.m @@ -34,6 +34,10 @@ % '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) % 'ece' = electron cyclotron emission % 'sxr' = soft x-ray emission % 'sxR' = soft x-ray emission with varargout{1} option (requires varargin{4}!) @@ -203,6 +207,18 @@ 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='vtor'; +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='vtorfit'; +end +if ~isempty(strmatch(data_type_eff_noext,[{'vpol'} {'v_pol'} {'VPOL'} {'V_POL'}],'exact')) + data_type_eff_noext='vpol'; +end +if ~isempty(strmatch(data_type_eff_noext,[{'vpolfit'} {'vpolft'} {'v_polfit'} {'VPOLfit'} {'V_POLfit'}],'exact')) + data_type_eff_noext='vpolfit'; +end % some defaults @@ -226,14 +242,15 @@ liuqeFBTE=[{'\results::i_p'} {'\results::z_axis'} {'\results::r_axis'} {'\result % 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'}]; + {'neft'} {'neft:trial'} {'teft:trial'} {'neftav:trial'} {'teftav:trial'} {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'} {'neint'} ... + {'vtor'} {'vtorfit'} {'vpol'} {'vpolfit'}]; % all keywords and corresponding case to run below TCVkeywrdall=[{'Ip'} {'B0'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'rhovol'} {'qrho'} {'q95'} {'kappa'} ... {'delta'} {'deltatop'} {'deltabot'} {'neint'} ... {'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'}]; + {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'} {'pgyro'} {'jtor'} {'vtor'} {'vtorfit'} {'vpol'} {'vpolfit'}]; TCVsig.iip=strmatch('Ip',TCVkeywrdall,'exact'); TCVsig.iB0=strmatch('B0',TCVkeywrdall,'exact'); TCVsig.izmag=strmatch('zmag',TCVkeywrdall,'exact'); @@ -277,6 +294,10 @@ 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.ivtor=strmatch('vtor',TCVkeywrdall,'exact'); +TCVsig.ivtorfit=strmatch('vtorfit',TCVkeywrdall,'exact'); +TCVsig.ivpol=strmatch('vpol',TCVkeywrdall,'exact'); +TCVsig.ivpolfit=strmatch('vpolfit',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 @@ -310,6 +331,10 @@ 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.ivtor)=TCVkeywrdall(TCVsig.ivtor); +TCVkeywrdcase(TCVsig.ivtorfit)=TCVkeywrdall(TCVsig.ivtorfit); +TCVkeywrdcase(TCVsig.ivpol)=TCVkeywrdall(TCVsig.ivpol); +TCVkeywrdcase(TCVsig.ivpolfit)=TCVkeywrdall(TCVsig.ivpolfit); % 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 @@ -465,18 +490,29 @@ switch TCVkeywrdcase{index} elseif length(tracetdi.dim{2})==41 ix=2; trace.t=tracetdi.dim{1}; - else - % time array unknow, 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 + 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 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 - disp(['assumes time array has more digits, so x from dim{' num2str(ix) '}']) end end if length(tracetdi.dim)==2 @@ -1170,6 +1206,102 @@ switch TCVkeywrdcase{index} mdsclose; error=0; + %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& + case {'vtor'} + mdsopen(shot); + nodenameeff=[{'\results::cxrs:vi_tor'} endstr]; + tracetdi=tdi(nodenameeff{:}); + trace.data=tracetdi.data; + trace.x=tracetdi.dim{1}; + trace.t=tracetdi.dim{2}; + 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 + nodenameeff=[{'\results::cxrs:vi_tor:err'} endstr]; + tracetdi=tdi(nodenameeff{:}); + trace.std = tracetdi.data; + trace.std_rho = tracetdi.dim{1}; + trace.std_t = tracetdi.dim{2}; + mdsclose; + error=0; + + %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& + case {'vtorfit'} + mdsopen(shot); + nodenameeff=[{'\results::cxrs:proffit:vi_tor'} endstr]; + tracetdi=tdi(nodenameeff{:}); + trace.data=tracetdi.data; + trace.x=tracetdi.dim{1}; + trace.t=tracetdi.dim{2}; + 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 + nodenameeff=[{'\results::cxrs:proffit:vi_tor:err'} endstr]; + tracetdi=tdi(nodenameeff{:}); + trace.std = tracetdi.data; + trace.std_rho = tracetdi.dim{1}; + trace.std_t = tracetdi.dim{2}; + mdsclose; + error=0; + + case {'vpol'} + mdsopen(shot); + nodenameeff=[{'\results::cxrs:vi_pol'} endstr]; + tracetdi=tdi(nodenameeff{:}); + trace.data=tracetdi.data; + trace.x=tracetdi.dim{1}; + trace.t=tracetdi.dim{2}; + 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 + nodenameeff=[{'\results::cxrs:vi_pol:err'} endstr]; + tracetdi=tdi(nodenameeff{:}); + trace.std = tracetdi.data; + trace.std_rho = tracetdi.dim{1}; + trace.std_t = tracetdi.dim{2}; + mdsclose; + error=0; + + %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& + case {'vpolfit'} + mdsopen(shot); + nodenameeff=[{'\results::cxrs:proffit:vi_pol'} endstr]; + tracetdi=tdi(nodenameeff{:}); + trace.data=tracetdi.data; + trace.x=tracetdi.dim{1}; + trace.t=tracetdi.dim{2}; + 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 + nodenameeff=[{'\results::cxrs:proffit:vi_pol:err'} endstr]; + tracetdi=tdi(nodenameeff{:}); + trace.std = tracetdi.data; + trace.std_rho = tracetdi.dim{1}; + trace.std_t = tracetdi.dim{2}; + mdsclose; + error=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']); diff --git a/crpptbx/gdat.m b/crpptbx/gdat.m index 242b82eaac33cc498960f221e84330fcf9cbd365..e0cee02e54154e548625f15d9c28837305cb5846 100644 --- a/crpptbx/gdat.m +++ b/crpptbx/gdat.m @@ -39,6 +39,10 @@ function [trace,error,varargout] = gdat(shot,data_type,varargin) % 'sxr' = soft x-ray emission % 'sxR' = soft x-ray emission with varargout{1} option (requires varargin{4}!) % 'MPX' = soft x-ray from wire chambers +% 'vtor' = toroidal rotation, raw data and arror bars +% 'vtorfit' = fitted toroidal rotation profiles (and error bars) +% 'vpol' = poloidal rotation, raw data and arror bars +% 'vpolfit' = fitted poloidal rotation profiles (and error bars) % % JET % Special case compatible with old gdat.m allows (JET related): gdat(51994,'ppf','efit/xip',1) @@ -70,6 +74,7 @@ function [trace,error,varargout] = gdat(shot,data_type,varargin) % trace.x: space of reference % trace.dim: cell array of grids, trace.dim{1}, {2}, ... % trace.dimunits: units of dimensions +% trace.std: if error bars are available % % error: error in loading signal (0=> OK, >=1 => error) % @@ -136,30 +141,38 @@ end % PLOT DATA (if required) if doplot==1 & length(trace.data)>1 & ~ischar(trace.data) - figure;zoom on - if length(size(trace.data))<=2 - plot(trace.t,trace.data); - ylabel(data_type) - else - for idim=1:length(trace.dim) - if length(trace.t)==length(trace.dim{idim}); idim_t=idim; end - end - if idim_t<=2 - plot(trace.t,trace.data(:,:,floor(end/2))); - ylabel([data_type '(:,:,floor(end/2))']) - elseif idim_t==3; - plot(trace.t,reshape(trace.data(:,floor(end/2),:),length(trace.dim{1}),length(trace.t))); - ylabel([data_type '(:,floor(end/2),:)']) + try + figure;zoom on + if length(size(trace.data))<=2 + plot(trace.t,trace.data); + ylabel(data_type) + else + for idim=1:length(trace.dim) + if length(trace.t)==length(trace.dim{idim}); idim_t=idim; end + end + if idim_t<=2 + plot(trace.t,trace.data(:,:,floor(end/2))); + ylabel([data_type '(:,:,floor(end/2))']) + elseif idim_t==3; + plot(trace.t,reshape(trace.data(:,floor(end/2),:),length(trace.dim{1}),length(trace.t))); + ylabel([data_type '(:,floor(end/2),:)']) + end end + xlabel('time [s]') + title([machine ' ' num2str(shot)]) + grid on + catch + disp(' error in plotting part, most probably because could not guess time dimension correctly. To check') end - xlabel('time [s]') - title([machine ' ' num2str(shot)]) - grid on elseif doplot==-1 - hold on - if length(size(trace.data))<=2 - plot(trace.t,trace.data,'r'); - else - plot(trace.t,trace.data(:,:,1),'--'); + try + hold on + if length(size(trace.data))<=2 + plot(trace.t,trace.data,'r'); + else + plot(trace.t,trace.data(:,:,1),'--'); + end + catch + disp(' error in plotting part, most probably because could not guess time dimension correctly. To check') end end