diff --git a/TCV/loadTCVdata.m b/TCV/loadTCVdata.m
index 5486788d3b540dbda659391985a9f30216fc4b14..c48b2446030071e57c5faa7d05f53afd02c2ae03 100644
--- a/TCV/loadTCVdata.m
+++ b/TCV/loadTCVdata.m
@@ -19,8 +19,10 @@
% '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
-% 'profnerho' = ne smoothed or fitted , vs (rho,t) (from Thomson fit)
-% 'profterho' = te smoothed or fitted , vs (rho,t) (from Thomson fit)
+% '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)
% '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)
@@ -47,6 +49,9 @@
% 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:
%
@@ -103,6 +108,9 @@ 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
@@ -157,7 +165,7 @@ liuqe23=[{'\results::i_p'} {'\results::z_axis'} {'\results::r_axis'} {'\results:
% all keywords and corresponding case to run below
TCVkeywrdall=[{'Ip'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'rhovol'} {'qrho'} {'q95'} {'kappa'} ...
{'delta'} {'deltatop'} {'deltabot'} {'neint'} ...
- {'ne'} {'te'} {'nerho'} {'terho'} {'profnerho'} {'profterho'} ...
+ {'ne'} {'te'} {'nerho'} {'terho'} {'nerhozshift'} {'terhozshift'} {'profnerho'} {'profterho'} ...
{'neft'} {'teft'} {'neftav'} {'teftav'} ...
{'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'}];
TCVsig.iip=strmatch('Ip',TCVkeywrdall,'exact');
@@ -178,6 +186,8 @@ TCVsig.ine=strmatch('ne',TCVkeywrdall,'exact');
TCVsig.ite=strmatch('te',TCVkeywrdall,'exact');
TCVsig.inerho=strmatch('nerho',TCVkeywrdall,'exact');
TCVsig.iterho=strmatch('terho',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');
@@ -202,6 +212,10 @@ TCVkeywrdcase(TCVsig.ine)=TCVkeywrdall(TCVsig.ine); % special as dimensions from
TCVkeywrdcase(TCVsig.ite)=TCVkeywrdall(TCVsig.ite); % idem
TCVkeywrdcase(TCVsig.inerho)=TCVkeywrdall(TCVsig.inerho); % idem
TCVkeywrdcase(TCVsig.iterho)=TCVkeywrdall(TCVsig.iterho); % 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);
@@ -228,8 +242,8 @@ 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.iprofnerho)={'\results::th_prof_ne'};
-TCVsiglocation(TCVsig.iprofterho)={'\results::th_prof_te'};
+%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;
@@ -329,17 +343,26 @@ switch TCVkeywrdcase{index}
trace.t=tracetdi.dim{max(1,TCVsigtimeindx(index))};
ix=3-TCVsigtimeindx(index); % works only for 2D arrays
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;
+ % if dim=41 assumes is x (usual for TCV)
+ if length(tracetdi.dim{1})==41
+ ix=1;
trace.t=tracetdi.dim{2};
- else
+ 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
+ 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
if length(tracetdi.dim)==2
trace.x=tracetdi.dim{ix};
@@ -368,7 +391,8 @@ switch TCVkeywrdcase{index}
trace.dim=tracetdi.dim;
trace.dimunits=tracetdi.dimunits;
end
- if isfield(tracetdi,'units');
+ % 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
@@ -399,7 +423,8 @@ switch TCVkeywrdcase{index}
trace.dimunits=[{'Z [m]'} ; {'time [s]'}];
trace.x=z;
trace.t=time;
- if isfield(tracetdi,'units');
+ % 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')
@@ -423,7 +448,6 @@ switch TCVkeywrdcase{index}
% add correct dimensions
time=mdsdata('\results::thomson:times');
% construct rho mesh
- liuqe_ext
psi_max=tdi(['\results::thomson:psi_max' liuqe_ext]);
psiscatvol=tdi(['\results::thomson:psiscatvol' liuqe_ext]);
for ir=1:length(psiscatvol.dim{2})
@@ -433,9 +457,104 @@ switch TCVkeywrdcase{index}
trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
trace.x=rho;
trace.t=time;
- if isfield(tracetdi,'units');
+ % 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 {'nerhozshift','terhozshift'}
+ % ne or Te from Thomson data on rho=sqrt(psi_normalised) mesh: (rho,t)
+ % allow for z shift of equil
+ if nargin>=2 & ~isempty(varargin{1})
+ zshift=varargin{1};
+ else
+ zshift=0.;
+ end
+ mdsopen(shot);
+ if strcmp(TCVkeywrdcase{index},'nerhozshift')
+ 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');
+ end
+ trace.data=tracetdi.data'; % Thomson data as (t,z)
+ trace.std=tracestd.data';
+ trace.name=[num2str(shot) ';' nodenameeff];
+ % add correct dimensions
+ time=mdsdata('\results::thomson:times');
+ % construct rho mesh
+ psi_max=tdi(['\results::thomson:psi_max' liuqe_ext]);
+ psiscatvol=tdi(['\results::thomson:psiscatvol' liuqe_ext]);
+ 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
+ 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
+ mdsclose
+
+ %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
+ case {'profnerho','profterho'}
+ % vol from psitbx
+ mdsopen(shot);
+ if strcmp(TCVkeywrdcase{index},'profnerho')
+ nodenameeff=['\results::THOMSON.PROFILES.AUTO:ne'];
+ end
+ if strcmp(TCVkeywrdcase{index},'profterho')
+ nodenameeff=['\results::THOMSON.PROFILES.AUTO:te'];
+ end
+ tracetdi=tdi(nodenameeff);
+ trace.data=tracetdi.data'; % error in dimensions for autofits
+ if ~isempty(tracetdi.dim)
+ trace.x=tracetdi.dim{1};
+ trace.t=tracetdi.dim{2};
+ else
+ 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
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
@@ -449,7 +568,8 @@ switch TCVkeywrdcase{index}
trace.t=tracetdi.dim{2};
trace.dim=[{trace.x};{trace.t}];
trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
- if isfield(tracetdi,'units');
+ % 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];
@@ -481,7 +601,8 @@ switch TCVkeywrdcase{index}
trace.dim=[{trace.x};{trace.t}];
end
trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
- if isfield(tracetdi,'units');
+ % 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 liuqe_ext];
@@ -510,7 +631,8 @@ switch TCVkeywrdcase{index}
trace.t=tracetdi.dim{2};
trace.dim=[{trace.x};{trace.t}];
trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
- if isfield(tracetdi,'units');
+ % 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 liuqe_ext];
@@ -596,7 +718,8 @@ switch TCVkeywrdcase{index}
radius.t=trace.t;
varargout{1}={radius};
error=0;
- if isfield(tracetdi,'units');
+ % 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
@@ -638,7 +761,8 @@ switch TCVkeywrdcase{index}
trace.t=tracetdi.dim{1};
trace.dim=tracetdi.dim;
trace.dimunits={'time [s]'};
- if isfield(tracetdi,'units');
+ % 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}];
@@ -656,7 +780,8 @@ switch TCVkeywrdcase{index}
trace.t=tracetdi.dim{1};
trace.dim=tracetdi.dim;
trace.dimunits={'time [s]'};
- if isfield(tracetdi,'units');
+ % 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}];