function [trace,error,varargout]=loadTCVdata(shot,data_type,varargin)
%
% list of data_type currently available:
% 'Ip' = current
% 'zmag' = vertical position of the center of the plasma (magnetic axis)
% 'rmag' = radial position of the center of the plasma
% 'rcont' = R of plama flux surfaces
% 'zcont' = Z of plama flux surfaces
% 'vol' = volume of flux surfaces
% 'qrho' = q profile on rho mesh
% 'neint' = line-integrated electron density [m/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'= ne profile on (rho=sqrt(psi),time) mesh.Note rho is a 2D array as depends on time. ADD error bars in .std
% 'terho'= 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 with .std=error bars, vs (rho,t) (from Thomson fit)
% 'profterho' = te smoothed or fitted with .std=error bars, vs (rho,t) (from Thomson fit)
% 'neft' = ne fitted from data on rho mesh (from proffit.local_time:neft)
% '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)
% 'teftav' = te fitted from averaged over time data on rho mesh (from proffit.avg_time:teft)
% '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
%
% 'xx_2 or xx_3' for Liuqe2 or 3: 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}: channel status 1= unread yet, 0= read
% (for traces with many channel, enables to load additional channels,
% like SXR, ECE, etc.)
% varargin{2}: [i1 i2] : if not empty, assumes need many chords from i1 to i2
% varargin{3}: zmag for varargout{1} computation
%
% 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
%
% 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]=loadTCVdata(shot,'Ip',1);
% [sxr,error,R]=loadTCVdata(shot,'sxR',1);
%
varargout{1}=cell(1,1);
error=1;
% To allow multiple ways of writing a specific keyword, use data_type_eff within this routine
data_type_eff=data_type;
i_23=0;
liuqe_ext='';
if strcmp(data_type_eff(end-1:end),'_2') | strcmp(data_type_eff(end-1:end),'_3')
i_23=2;
liuqe_ext=data_type_eff(end-1:end);
end
% use keyword without eventual _2 or _3 extension to check for multiple possibilities
data_type_eff_noext=data_type_eff(1:end-i_23);
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,[{'Te'} {'t_e'} {'TE'} {'T_e'}],'exact'))
data_type_eff_noext='te';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Ne'} {'n_e'} {'NE'} {'N_e'}],'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,[{'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,[{'Rmag'}],'exact'))
data_type_eff_noext='rmag';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Zmag'}],'exact'))
data_type_eff_noext='zmag';
end
% some defaults
% nodes which have _2 and _3 equivalence, related to simpletdi case
liuqe23=[{'\results::i_p'} {'\results::z_axis'} {'\results::r_axis'} {'\results::q_psi'} ...
{'\results::beta_tor'} {'\results::q_95'} {'\results::l_i'} {'\results::delta_95'} ...
{'\results::kappa_95'} {'\results::r_contour'} {'\results::z_contour'} {'\results::psi_axis'} ...
{'\results::thomson:psiscatvol'} {'\results::thomson:psi_max'} {'\results::psitbx:vol'}];
% all keywords and corresponding case to run below
TCVkeywrdall=[{'Ip'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'qrho'} {'neint'} ...
{'ne'} {'te'} {'nerho'} {'terho'} {'profnerho'} {'profterho'} ...
{'neft'} {'teft'} {'neftav'} {'teftav'} ...
{'sxr'} {'sxR'} {'ece'}];
TCVsig.iip=strmatch('Ip',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.iqrho=strmatch('qrho',TCVkeywrdall,'exact');
TCVsig.ineint=strmatch('neint',TCVkeywrdall,'exact');
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.iprofnerho=strmatch('profnerho',TCVkeywrdall,'exact');
TCVsig.iprofterho=strmatch('profterho',TCVkeywrdall,'exact');
TCVsig.ineft=strmatch('neft',TCVkeywrdall,'exact');
TCVsig.iteft=strmatch('teft',TCVkeywrdall,'exact');
TCVsig.ineftav=strmatch('neftav',TCVkeywrdall,'exact');
TCVsig.iteftav=strmatch('teftav',TCVkeywrdall,'exact');
TCVsig.isxr=strmatch('sxr',TCVkeywrdall,'exact');
TCVsig.isxR=strmatch('sxR',TCVkeywrdall,'exact');
TCVsig.iece=strmatch('ece',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.iqrho)=TCVkeywrdall(TCVsig.iqrho); % special as liuqe q_psi on psi
TCVkeywrdcase(TCVsig.ine)=TCVkeywrdall(TCVsig.ine); % special as dimensions from other nodes
TCVkeywrdcase(TCVsig.ite)=TCVkeywrdall(TCVsig.ite); % idem
TCVkeywrdcase(TCVsig.inerho)=TCVkeywrdall(TCVsig.inerho); % idem
TCVkeywrdcase(TCVsig.iterho)=TCVkeywrdall(TCVsig.iterho); % idem
TCVkeywrdcase(TCVsig.iprofnerho)=TCVkeywrdall(TCVsig.iprofnerho); % special as add .std
TCVkeywrdcase(TCVsig.iprofterho)=TCVkeywrdall(TCVsig.iprofterho); % idem
TCVkeywrdcase(TCVsig.isxr)=TCVkeywrdall(TCVsig.isxr);
TCVkeywrdcase(TCVsig.isxR)=TCVkeywrdall(TCVsig.isxR);
TCVkeywrdcase(TCVsig.iece)=TCVkeywrdall(TCVsig.iece);
% 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.ivol)={'\results::psitbx:vol'}; TCVsigtimeindx(TCVsig.ivol)=2;
TCVsiglocation(TCVsig.ineint)={'\results::fir:lin_int_dens'};
TCVsiglocation(TCVsig.ineft)={'\results::proffit.local_time:neft'}; TCVsigtimeindx(TCVsig.ineft)=2;
TCVsiglocation(TCVsig.iteft)={'\results::proffit.local_time:teft'}; TCVsigtimeindx(TCVsig.iteft)=2;
TCVsiglocation(TCVsig.ineftav)={'\results::proffit.avg_time:neft'}; TCVsigtimeindx(TCVsig.ineftav)=2;
TCVsiglocation(TCVsig.iteftav)={'\results::proffit.avg_time:teft'}; TCVsigtimeindx(TCVsig.iteftav)=2;
% 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 isempty(index)
disp('********************')
disp('trace not yet registered.')
disp('If standard data, ask andrea.scarabosio@epfl.ch or 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('********************')
% 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(1:end-i_23),TCVkeywrdall,'exact');
if isempty(index)
disp(' ')
disp('********************')
disp(['no such keyword: ' data_type_eff(1:end-i_23)])
disp(' ')
disp('Available keywords:')
TCVkeywrdall(:)
disp('********************')
trace.data=[];
trace.t=[];
trace.x=[];
return
end
end
disp(['loading' ' ' data_type_eff ' from TCV shot #' num2str(shot)]);
disp(['case ' TCVkeywrdcase{index}])
if i_23==2 & isempty(strmatch(TCVsiglocation(index),liuqe23,'exact'))
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;
i1=1;
i2=20;
if nargineff>=3
for i=1:length(varargin)
if ~isempty(varargin{i})
if isstruct(varargin{i})
zmag=varargin{i};
elseif size(varargin{i},2)>2
status=varargin{i};
else
i1 =varargin{i}(1);
i2 =varargin{i}(2);
end
end
end
end
switch TCVkeywrdcase{index}
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case 'simpletdi'
% load TCV other data
mdsopen(shot);
nodenameeff=[TCVsiglocation{index} liuqe_ext];
% test if node exists
error=1;
if eval(['~mdsdata(''node_exists("\' nodenameeff '")'')'])
disp(['node ' nodenameeff ' does not exist for shot = ' num2str(shot)])
trace.data=[];
trace.x=[];
trace.t=[];
trace.dim=[];
trace.dimunits=[];
return
elseif eval(['mdsdata(''getnci("\' nodenameeff ':foo","length")'')==0'])
disp(['no data for node ' nodenameeff ' for shot = ' num2str(shot)])
trace.data=[];
trace.x=[];
trace.t=[];
trace.dim=[];
trace.dimunits=[];
return
end
tracetdi=tdi(nodenameeff);
mdsclose(shot)
if isempty(tracetdi.data) | isnan(tracetdi.data)
disp(['node ' nodenameeff ' is empty for shot = ' num2str(shot)])
trace.data=[];
trace.x=[];
trace.t=[];
trace.dim=[];
trace.dimunits=[];
return
end
trace.data=tracetdi.data;
trace.x=[]; % so that appears before trace.t in list of sub-structures
if TCVsigtimeindx(index)>0
trace.t=tracetdi.dim{TCVsigtimeindx(index)};
ix=3-TCVsigtimeindx(index); % works only for 2D arrays
else
trace.t=tracetdi.dim{1};
% 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; 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};
% 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
trace.x=[];
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(1)=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;
trace.x=[];
trace.t=[];
end
error=0;
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case {'ne','te'}
% ne or Te from Thomson data on raw z mesh vs (z,t)
mdsopen(shot);
if strcmp(TCVkeywrdcase{index},'ne')
tracetdi=tdi('\results::thomson:ne');
tracestd=tdi('\results::thomson:ne:error_bar');
else
tracetdi=tdi('\results::thomson:te');
tracestd=tdi('\results::thomson:te:error_bar');
end
trace.data=tracetdi.data'; % Thomson data as (t,z)
trace.std=tracestd.data;
% add correct dimensions
time=mdsdata('\results::thomson:times');
z=mdsdata('\diagz::thomson_set_up:vertical_pos');
trace.dim=[{z};{time}];
trace.dimunits=[{'Z [m]'} ; {'time [s]'}];
trace.x=z;
trace.t=time;
mdsclose
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case {'nerho','terho'}
% ne or Te from Thomson data on rho=sqrt(psi_normalised) mesh: (rho,t)
mdsopen(shot);
if strcmp(TCVkeywrdcase{index},'nerho')
tracetdi=tdi('\results::thomson:ne');
tracestd=tdi('\results::thomson:ne:error_bar');
else
tracetdi=tdi('\results::thomson:te');
tracestd=tdi('\results::thomson:te:error_bar');
end
trace.data=tracetdi.data'; % Thomson data as (t,z)
trace.std=tracestd.data;
% 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]);
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;
mdsclose
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case {'qrho'}
% q profile on psi from liuqe
mdsopen(shot);
tracetdi=tdi(['\results::q_psi' liuqe_ext]);
trace.data=tracetdi.data;
trace.x=sqrt(tracetdi.dim{1}/10);
trace.t=tracetdi.dim{2};
trace.dim=[{trace.x};{trace.t}];
trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
mdsclose
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case {TCVkeywrdcase{TCVsig.iprofnerho},TCVkeywrdcase{TCVsig.iprofnerho}}
% Thomson profiles with error bars
mdsopen(shot);
if i_23==2
disp('cannot have profile on Liuqe>1, ask roland.behn@epfl.ch')
return
end
if strcmp(TCVkeywrdcase{index},TCVkeywrdcase{TCVsig.iprofnerho})
tracetdi=tdi('\results::th_prof_ne');
tracestd=tdi('\results::thomson:ne:error_bar');
else
tracetdi=tdi('\results::th_prof_te');
tracestd=tdi('\results::thomson:te:error_bar');
end
trace.data=tracetdi.data'; % Thomson data transposed
trace.x=tracetdi.dim{2};
trace.t=tracetdi.dim{1};
trace.dim=[{trace.x};{trace.t}];
trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
trace.std=tracestd.data;
mdsclose
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case {'sxr','sxR'}
% load TCV soft x-ray data
% 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 isempty(zmag)
zmag=loadTCVdata(shot,['zmag' liuqe_ext]);
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.t=t';
radius.data={VsxrTCVradius(data,xchord,ychord)};
varargout{1}={radius};
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
for i=1:(20*icamera*ones(10,1))
trace.t(:,i)=t';
end
trace.data=xtomo_signal';
trace.dim{1}={trace.t};
trace.dimunits={'time [s]'};
error=0;
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case 'ece'
% load TCV ECE data
% Status=1 => Not Read Yet
if ~isempty(find(status == 1))
[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));
for i=1:size(TE_ECE,2)
trace.t(:,i)=T(a:b);
end
trace.data=TE_ECE(a:b,:);
radius.t=trace.t;
radius.data=R(a:b,:);
trace.dim{1}={trace.t};
trace.dimunits={'time [s]'};
varargout{1}={radius};
error=0;
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
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);
signal=get_mds_mio('MPX',[t_1 t_2]);
mdsclose(shot)
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]'};
[xchord,ychord]=mpx_geometry;
varargout{1}={VsxrTCVradius(zmag.data,xchord,ychord)};
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']);
end