Newer
Older
function [trace,error,varargout]=loadAUGdata(shot,data_type,varargin)
%
% data_type:
% 'Ip' = current
% 'zmag' = vertical position of the center of the plasma (magnetic axis)
% 'rmag' = radial position of the center of the plasma
% 'sxr' = soft x-ray emission
% 'sxR' = soft x-ray emission with varargout{1} option (requires varargin{5}!)
% 'SXB' = soft x-ray emission from (by default camera J) SXB/J_xx camera (sxb, sxB, etc all work)
% 'SXF' = soft x-ray emission from (by default camera J) SXF/I_xx camera (sxf, sxF, etc all work)
% 'SSX_H' = from SSX, H channel, same for G, I, etc
% 'SSX' : gets the SSX_G by default at this stage
% gdat(15133,'MAG/Ipa',1,'AUG') % uses experiment=AUGD shotfiles per default
% gdat(15133,'MAG/Ipa',1) (sufficient at AUG since 'AUG' per defaut, same with gdat(15133,'ip',1)
% gdat(15133,'AUGD/MAG/Ipa',1,'AUG') % to specify experiment explicitely like in:
% gdat(30230,'ECED/RMD/Trad-A',1,'AUG') %
%
% INPUT:
% shot: shot number
% data_type: type of the required data: 'diag_name/sig_name'
%
% examples:
% data_type='SXR/B', 'TOT/beta_N', 'SXB/J_053', 'SXB'
% data_type='POT/ELMa-Han', 'MOD/OddNAmp', 'MOD/EvenNAmp', 'TOT/PNBI_TOT', 'TOT/P_TOT'
% 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
% varargin{4}: time range [t1 t2] (to limit data collected)
% varargin{5}: nth to keep only nth time points of traces
% varargin{6}: camera to use: 'B' (default), 'A',
%
% OUTPUT:
% trace.data: data structure
% trace.t: time of reference
% trace.x: space of reference
% .... others related to data
% error: error in loading signal (0=> OK, 1=> error)
%
% Additional Output arguments depending on data_type
%
% data_type=sxR:
% varargout{1}: intersection of the view lines with magnetic axis
%
% functions needed: SF2ML or mdsplus routines
%
% Example:
% [ip,error]=loadAUGdata(shot,'ip');
% [ip,error]=loadAUGdata(shot,'MAG/Ipi');
% [n2,error]=loadAUGdata(shot,'MOD/EvenNAmp');
%
% List of keywords (can be used in adition to 'DIAG/tracenam'), with comments when not obvious:
%
% Ip :
% b0 :
% zmag :
% rmag :
% rgeo :
% zgeo :
% q0 :
% q95 :
% kappa :
% delta :
% deltatop :
% deltabot :
% neint :
% neterho :
% cxrs : CXRS structure from CEZ with vrot, Ti, ...
% cxrs_rho : same as cxrs but project on rho as well (using 'equil' information)
% equil : equilibrium structure, rhopol, rhotor, rhovol, etc from EQI
% equil_fpp : as equil but from FPP
% equil_eqm : as equil but from EQM
% equil_eqr : as equil but from EQR
% sxr : from SXR/A or B (from old stuff, not sure still OK)
% sxR : from SXR/A or B adding R of chords (from old stuff, not sure still OK)
% sxb : 'SXB/J' chords
% ece :
% eced :
% eced_rho :
% Halpha :
% pgyro : for each gyrotrons, power, freq, etc (ask for more)
% powers : power traces for each sources
%
shotfile_exp = 'AUGD';
% To allow multiple ways of writing a specific keyword, use data_type_eff within this routine
data_type_eff=data_type;
if size(data_type,1)==1
% there might be "/" in the name, add backslash to protect it, thus remove in the name afterwards:
inotok=findstr('\/',data_type);
iok=regexp(data_type,'[^\\]/');
iok=iok+1;
if ~isempty(inotok)
% remove \/ and construct index of real separators "/"
for ij=1:length(inotok)
ijk=find(iok>inotok(ij));
if ~isempty(ijk)
iok(ijk) = iok(ijk) - 1;
end
end
end
inobackslash=regexp(data_type,'[^\\]');
data_type=data_type(inobackslash);
index_slash = iok;
% i=findstr('/',data_type);
i = index_slash;
if length(i)==1
data_type_eff=[{data_type(1:i(1)-1)} ; {data_type(i(1)+1:end)}];
elseif length(i)==2
% assumes given a la 'AUGD/MAG/Ipi' or 'ECED/RMP/Trad-A'
data_type_eff=[{data_type(i(1)+1:i(2)-1)} ; {data_type(i(2)+1:end)}];
shotfile_exp = data_type(1:i(1)-1);
disp(['more / than expected in tracename: length(i)= ' num2str(length(i))])
data_type
Olivier Sauter
committed
elseif isempty(data_type)
data_type_eff = ' ';
end
i_efitm=0;
i_ext=length(data_type_eff)+1;
name_ext='';
if size(data_type_eff,1)==1
data_type_eff_noext=data_type_eff(1:i_ext-1);
if ~isempty(strmatch(data_type_eff_noext,[{'ip'} {'i_p'} {'IP'} {'iP'} {'xip'}],'exact'))
if ~isempty(strmatch(data_type_eff_noext,[{'B0'}],'exact'))
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,[{'Terho'}],'exact'))
data_type_eff_noext='terho';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'neterho'}],'exact'))
data_type_eff_noext='neterho';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'cxrs'} {'vrot'} {'ti'}],'exact'))
data_type_eff_noext='cxrs'; % load full CEZ structure
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'cxrs_rhos'} {'cxrs_rho'} {'cxrsrho'} {'cxrsrhos'}],'exact'))
data_type_eff_noext='cxrs_rho'; % load full CEZ structure
if ~isempty(strmatch(data_type_eff_noext,[{'SXR'}],'exact'))
data_type_eff_noext='sxr';
end
if ~isempty(strmatch(data_type_eff_noext,[{'SXB'} {'sxb'} {'Sxb'} {'sXb'} {'sxB'} {'SXb'}],'exact'))
data_type_eff_noext='sxb';
end
if ~isempty(strmatch(data_type_eff_noext,[{'SXF'} {'sxf'} {'Sxf'} {'sXf'} {'sxf'} {'SXf'}],'exact'))
data_type_eff_noext='sxf';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'ssx_g'} {'ssx'}],'exact'))
data_type_eff_noext='ssx_g';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'ssx_h'} {'ssx'}],'exact'))
data_type_eff_noext='ssx_h';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'ssx_i'} {'ssx'}],'exact'))
data_type_eff_noext='ssx_i';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'ssx_j'} {'ssx'}],'exact'))
data_type_eff_noext='ssx_j';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'ssx'}],'exact'))
data_type_eff_noext='ssx';
end
if ~isempty(strmatch(data_type_eff_noext,[{'ECE'}],'exact'))
data_type_eff_noext='ece';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'ece_rho'}],'exact'))
data_type_eff_noext='ece_rho';
end
if ~isempty(strmatch(upper(data_type_eff_noext),[{'ECED'}],'exact'))
data_type_eff_noext='eced';
end
if ~isempty(strmatch(upper(data_type_eff_noext),[{'ECED_RHO'}],'exact'))
data_type_eff_noext='eced_rho';
end
if ~isempty(strmatch(upper(data_type_eff_noext),[{'ECED_RMD'}],'exact'))
data_type_eff_noext='eced_rmd';
end
if ~isempty(strmatch(data_type_eff_noext,[{'VOL'} {'volume'}],'exact'))
data_type_eff_noext='vol';
end
if ~isempty(strmatch(data_type_eff_noext,[{'q_0'} {'Q0'}],'exact'))
data_type_eff_noext='q0';
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(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,[{'qrho_FPP'}],'exact'))
data_type_eff_noext='qrho_fpp';
end
if ~isempty(strmatch(data_type_eff_noext,[{'equil_FPP'}],'exact'))
data_type_eff_noext='equil_fpp';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'equil_eqm'}],'exact'))
data_type_eff_noext='equil_eqm';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'equil_eqr'}],'exact'))
data_type_eff_noext='equil_eqr';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Rmag'}],'exact'))
data_type_eff_noext='rmag';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Rgeo'}],'exact'))
data_type_eff_noext='rgeo';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Zmag'}],'exact'))
data_type_eff_noext='zmag';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Zgeo'}],'exact'))
data_type_eff_noext='zgeo';
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,[{'Ha'} {'ha'} {'Halpha'} {'halpha'}],'exact'))
data_type_eff_noext='Halpha';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'pgyro'} {'pec'} {'pech'} {'pecrh'} {'p_ec'} {'p_gyro'}],'exact'))
data_type_eff_noext='pgyro';
end
if ~isempty(strmatch(lower(data_type_eff_noext),[{'powers'} {'ptot'} {'ptots'}],'exact'))
data_type_eff_noext='powers';
end
else
i_ext=length(data_type_eff{2})+1;
name_ext='';
data_type_eff_noext=data_type_eff{2}(1:i_ext-1);
end
% all keywords and corresponding case to run below
AUGkeywrdall=[{'Ip'} {'b0'} {'zmag'} {'rmag'} {'rgeo'} {'zgeo'} {'rcont'} {'zcont'} {'vol'} {'qrho'} {'qrho_fpp'} {'q0'} {'q95'} {'kappa'} ...
{'delta'} {'deltatop'} {'deltabot'} {'neint'} {'ne'} {'te'} ...
{'nerho'} {'neterho'} {'terho'} {'cxrs'} {'cxrs_rho'} {'equil'} {'equil_fpp'} {'equil_eqm'} ...
{'equil_eqr'} {'sxr'} {'sxR'} {'sxb'} {'sxf'} {'ssx_g'} {'ssx_h'} {'ssx_i'} {'ssx_j'} {'ssx'} ...
{'ece'} {'ece_rho'} {'eced'} {'eced_rho'} {'eced_rmd'} {'Halpha'} {'pgyro'} {'powers'}];
AUGsig.ib0=strmatch('b0',AUGkeywrdall,'exact');
AUGsig.izmag=strmatch('zmag',AUGkeywrdall,'exact');
AUGsig.irmag=strmatch('rmag',AUGkeywrdall,'exact');
AUGsig.irgeo=strmatch('rgeo',AUGkeywrdall,'exact');
AUGsig.izgeo=strmatch('zgeo',AUGkeywrdall,'exact');
AUGsig.ircont=strmatch('rcont',AUGkeywrdall,'exact');
AUGsig.izcont=strmatch('zcont',AUGkeywrdall,'exact');
AUGsig.ivol=strmatch('vol',AUGkeywrdall,'exact');
AUGsig.iqrho=strmatch('qrho',AUGkeywrdall,'exact');
AUGsig.iqrho_fpp=strmatch('qrho_fpp',AUGkeywrdall,'exact');
AUGsig.iequil=strmatch('equil',AUGkeywrdall,'exact');
AUGsig.iequil_fpp=strmatch('equil_fpp',AUGkeywrdall,'exact');
AUGsig.iequil_eqm=strmatch('equil_eqm',AUGkeywrdall,'exact');
AUGsig.iequil_eqr=strmatch('equil_eqr',AUGkeywrdall,'exact');
AUGsig.iq0=strmatch('q0',AUGkeywrdall,'exact');
AUGsig.iq95=strmatch('q95',AUGkeywrdall,'exact');
AUGsig.ikappa=strmatch('kappa',AUGkeywrdall,'exact');
AUGsig.idelta=strmatch('delta',AUGkeywrdall,'exact');
AUGsig.ideltatop=strmatch('deltatop',AUGkeywrdall,'exact');
AUGsig.ideltabot=strmatch('deltabot',AUGkeywrdall,'exact');
AUGsig.ineint=strmatch('neint',AUGkeywrdall,'exact');
AUGsig.ine=strmatch('ne',AUGkeywrdall,'exact');
AUGsig.ite=strmatch('te',AUGkeywrdall,'exact');
AUGsig.inerho=strmatch('nerho',AUGkeywrdall,'exact');
AUGsig.iterho=strmatch('terho',AUGkeywrdall,'exact');
AUGsig.ineterho=strmatch('neterho',AUGkeywrdall,'exact');
AUGsig.icxrs=strmatch('cxrs',AUGkeywrdall,'exact');
AUGsig.icxrs_rho=strmatch('cxrs_rho',AUGkeywrdall,'exact');
AUGsig.isxr=strmatch('sxr',AUGkeywrdall,'exact');
AUGsig.isxR=strmatch('sxR',AUGkeywrdall,'exact');
AUGsig.isxb=strmatch('sxb',AUGkeywrdall,'exact');
AUGsig.isxf=strmatch('sxf',AUGkeywrdall,'exact');
AUGsig.issx_g=strmatch('ssx_g',AUGkeywrdall,'exact');
AUGsig.issx_h=strmatch('ssx_h',AUGkeywrdall,'exact');
AUGsig.issx_i=strmatch('ssx_i',AUGkeywrdall,'exact');
AUGsig.issx_j=strmatch('ssx_j',AUGkeywrdall,'exact');
AUGsig.issx=strmatch('ssx',AUGkeywrdall,'exact');
AUGsig.ieced=strmatch('eced',AUGkeywrdall,'exact');
AUGsig.iece_rho=strmatch('ece_rho',AUGkeywrdall,'exact');
AUGsig.ieced_rho=strmatch('eced_rho',AUGkeywrdall,'exact');
AUGsig.ieced_rmd=strmatch('eced_rmd',AUGkeywrdall,'exact');
AUGsig.iHalpha=strmatch('Halpha',AUGkeywrdall,'exact');
AUGsig.ipgyro=strmatch('pgyro',AUGkeywrdall,'exact');
AUGsig.ipowers=strmatch('powers',AUGkeywrdall,'exact');
% For each keyword, specify which case to use. As most common is 'simplereaddata', fill in with this and change
% only indices needed. Usually use name of case same as keyword name
AUGkeywrdcase=cell(size(AUGkeywrdall));
AUGkeywrdcase(:)={'simplereaddata'};
AUGkeywrdcase(AUGsig.iqrho)=AUGkeywrdall(AUGsig.iqrho); % special as efit q on psi
AUGkeywrdcase(AUGsig.iqrho_fpp)=AUGkeywrdall(AUGsig.iqrho_fpp); % special as efit q on psi
AUGkeywrdcase(AUGsig.iequil)=AUGkeywrdall(AUGsig.iequil); % special as efit q on psi
AUGkeywrdcase(AUGsig.iequil_fpp)=AUGkeywrdall(AUGsig.iequil_fpp); % special as efit q on psi
AUGkeywrdcase(AUGsig.iequil_eqm)=AUGkeywrdall(AUGsig.iequil_eqm); % special as efit q on psi
AUGkeywrdcase(AUGsig.iequil_eqr)=AUGkeywrdall(AUGsig.iequil_eqr); % special as efit q on psi
%AUGkeywrdcase(AUGsig.idelta)=AUGkeywrdall(AUGsig.idelta); % special as average of triu and tril
AUGkeywrdcase(AUGsig.ine)=AUGkeywrdall(AUGsig.ine); % special as adds error bars
AUGkeywrdcase(AUGsig.ite)=AUGkeywrdall(AUGsig.ite); % idem
AUGkeywrdcase(AUGsig.inerho)=AUGkeywrdall(AUGsig.inerho); % idem
AUGkeywrdcase(AUGsig.ineterho)=AUGkeywrdall(AUGsig.ineterho); % idem
AUGkeywrdcase(AUGsig.iterho)=AUGkeywrdall(AUGsig.iterho); % idem
AUGkeywrdcase(AUGsig.isxr)=AUGkeywrdall(AUGsig.isxr);
AUGkeywrdcase(AUGsig.isxR)=AUGkeywrdall(AUGsig.isxR);
AUGkeywrdcase(AUGsig.isxb)=AUGkeywrdall(AUGsig.isxb);
AUGkeywrdcase(AUGsig.isxf)=AUGkeywrdall(AUGsig.isxf);
AUGkeywrdcase(AUGsig.issx_g)=AUGkeywrdall(AUGsig.issx_g);
AUGkeywrdcase(AUGsig.issx_h)=AUGkeywrdall(AUGsig.issx_h);
AUGkeywrdcase(AUGsig.issx_i)=AUGkeywrdall(AUGsig.issx_i);
AUGkeywrdcase(AUGsig.issx_j)=AUGkeywrdall(AUGsig.issx_j);
AUGkeywrdcase(AUGsig.issx)=AUGkeywrdall(AUGsig.issx);
AUGkeywrdcase(AUGsig.iece)=AUGkeywrdall(AUGsig.iece);
AUGkeywrdcase(AUGsig.ieced)=AUGkeywrdall(AUGsig.ieced);
AUGkeywrdcase(AUGsig.iece_rho)=AUGkeywrdall(AUGsig.iece_rho);
AUGkeywrdcase(AUGsig.ieced_rho)=AUGkeywrdall(AUGsig.ieced_rho);
AUGkeywrdcase(AUGsig.icxrs)=AUGkeywrdall(AUGsig.icxrs);
AUGkeywrdcase(AUGsig.icxrs_rho)=AUGkeywrdall(AUGsig.icxrs_rho);
AUGkeywrdcase(AUGsig.ipgyro)=AUGkeywrdall(AUGsig.ipgyro); % idem
AUGkeywrdcase(AUGsig.ipowers)=AUGkeywrdall(AUGsig.ipowers); % idem
% 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
AUGsigtimeindx=ones(size(AUGkeywrdall));
% For the 'simplereaddata' cases, we need the full node in case gdat was called with full location directly
% for the other cases, leave this location empty
AUGsiglocation=cell(2,size(AUGkeywrdall,2));
AUGsiglocation(:)={''};
AUGsiglocation(:,AUGsig.iip)={'MAG'; 'Ipa'};
AUGsiglocation(:,AUGsig.ib0)={'FPC'; 'BTF'}; % at 1.65m (?)
AUGsiglocation(:,AUGsig.izmag)={'FPG'; 'Zmag'};
AUGsiglocation(:,AUGsig.irmag)={'FPG'; 'Rmag'};
AUGsiglocation(:,AUGsig.irgeo)={'FPG'; 'Rgeo'};
AUGsiglocation(:,AUGsig.izgeo)={'FPG'; 'Zgeo'};
AUGsiglocation(:,AUGsig.ircont)={'' ; ''}; AUGsigtimeindx(AUGsig.ircont)=2;
AUGsiglocation(:,AUGsig.izcont)={'' ; ''}; AUGsigtimeindx(AUGsig.izcont)=2;
AUGsiglocation(:,AUGsig.ivol)={''; ''};
AUGsiglocation(:,AUGsig.iq0)={'FPG'; 'q0'};
AUGsiglocation(:,AUGsig.iq95)={'FPG'; 'q95'};
AUGsiglocation(:,AUGsig.ikappa)={''; ''};
AUGsiglocation(:,AUGsig.ideltatop)={''; ''};
AUGsiglocation(:,AUGsig.ideltabot)={''; ''};
AUGsiglocation(:,AUGsig.ineint)={'DCN'; 'H-1'};
AUGsiglocation(:,AUGsig.iece)={'CEC'; 'Trad-A'};
AUGsiglocation(:,AUGsig.ieced)={'CEC'; 'Trad-A'}; % ECED
AUGsiglocation(:,AUGsig.iece_rho)={'CEC'; 'Trad-A'};
AUGsiglocation(:,AUGsig.ieced_rho)={'CEC'; 'Trad-A'}; % ECED
AUGsiglocation(:,AUGsig.ieced_rmd)={'RMD'; 'Trad-A'}; % ECED_RMD
AUGsiglocation(:,AUGsig.iHalpha)={'POT'; 'ELMa-Han'};
% For the 'simplereaddata' cases, we need the full node in case gdat was called with full location directly
% for the other cases, leave this location empty
AUGexplocation=cell(size(AUGkeywrdall,2),1);
AUGexplocation(:)={'AUGD'};
% cases with a "private" shotfile:
AUGexplocation(AUGsig.ieced)={'ECED'};
AUGexplocation(AUGsig.ieced_rho)={'ECED'};
Olivier Sauter
committed
if shot==-9
clear trace
for i=1:length(AUGkeywrdall)
fieldname_eff = lower(AUGkeywrdall{i});
keyword_eff = AUGkeywrdall{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 = 'MAG/Ipa';
% trace.data=[];
return
end
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
% initialize order of substructures and allows just a "return" if data empty
trace.data=[];
trace.x=[];
trace.t=[];
trace.dim=[];
trace.dimunits=[];
trace.name=[];
% find index of signal called upon
if size(data_type_eff,1)==2
% in case node name was given in 2 parts directly (as old way)
ii1=strmatch(data_type_eff(1),AUGsiglocation(1,:),'exact');
iiindex=strmatch(data_type_eff_noext,AUGsiglocation(2,ii1),'exact');
if ~isempty(iiindex)
index=ii1(iiindex);
else
index=[];
end
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{1} ' ' data_type_eff{2} ' ' num2str(shot) ' ' ...
% getenv('USER') ' AUG'' olivier.sauter@epfl.ch < /dev/null'])
disp('********************')
% temporarily add entry in arrays, so can work below
index=length(AUGkeywrdall)+1;
AUGkeywrdall(end+1)={'new'};
AUGkeywrdcase(end+1)={'simplereaddata'};
AUGsiglocation(1:2,end+1)=[data_type_eff(1) ; {data_type_eff_noext}];
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
AUGsigtimeindx(end+1)=0;
elseif ~strcmp(AUGkeywrdcase{index},'simplereaddata')
msgbox(['Problem in loadAUGdata with data_type_eff = ' char(data_type_eff(end)) ...
'. Full paths of nodes should only be for case simplereaddata'],'in loadAUGdata','error')
error('in loadAUGdata')
end
else
index=strmatch(data_type_eff_noext,AUGkeywrdall,'exact');
if isempty(index)
disp(' ')
disp('********************')
if iscell(data_type_eff)
disp(['no such keyword: ' data_type_eff{1} '/' data_type_eff{2}])
else
disp(['no such keyword: ' data_type_eff])
end
disp(' ')
disp('Available keywords:')
AUGkeywrdall(:)
disp('********************')
return
end
end
disp(' ')
if iscell(data_type_eff)
disp(['loading' ' ' data_type_eff{1} '/' data_type_eff{2} ' from AUG shot #' num2str(shot)]);
else
disp(['loading' ' ' data_type_eff ' from AUG shot #' num2str(shot)]);
end
disp(['case ' AUGkeywrdcase{index}])
disp(' ')
switch AUGkeywrdcase{index}
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case 'simplereaddata'
ppftype=AUGsiglocation{1,index};
if i_efitm;
tracename=['eftm' AUGsiglocation{2,index}(5:end) name_ext];
else
tracename=[AUGsiglocation{2,index} name_ext];
end
ij=find(tracename~='''');
[a,e]=rdaAUG_eff(shot,ppftype,tracename,shotfile_exp_eff);
% switch tracename
% special cases if traces do not exist for some shot or other
% end
trace=a;
clear error
error=e;
if isempty(trace.data)
trace.dimunits=[];
elseif length(size(trace.data))==1 | (length(size(trace.data))==2 & size(trace.data,2)==1)
trace.dim=[{trace.t}];
trace.dimunits={'time [s]'};
elseif length(size(trace.data))==2
trace.dim=[{trace.x} ; {trace.t}];
trace.dimunits=[{'R [m] or rho=sqrt(psi_norm)'} ; {'time [s]'}];
else
% disp('how to deal with 3D arrays?')
% find time dimension
ilentime=length(trace.t);
ij=find(size(trace.value)==ilentime);
if ij==1;
% as expected
trace.x = [1:size(trace.value,2)];
trace.dim=[{trace.t} ; {trace.x} ; {[1:size(trace.value,3)]}];
trace.dimunits=[{'time [s]'} ; {''} ; {''}];
else
trace.dim=[{[]} ; {[]} ; {[]}];
trace.dim{ij} = trace.t;
trace.dimunits=[{[]} ; {[]} ; {[]}];
trace.dimunits{ij}='time [s]';
end
end
trace.name=[num2str(shot) '/' ppftype '/' tracename];
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case {'sxr','sxR'}
% LOAD MULTI CHANNEL DATA
% load AUG soft x-ray data
if nargin>=3 & ~isempty(varargin{1})
starti=varargin{1}(1);
endi=varargin{1}(2);
else
starti=1;
endi=30;
end
if nargin>=4 & ~isempty(varargin{2})
status=varargin{2};
else
status=ones(endi-starti+1,1);
end
if nargin>=6 & ~isempty(varargin{4})
timerange=varargin{4};
else
if nargin>=7 & ~isempty(varargin{5})
nth=varargin{5};
else
nth=13;
end
if nargin>=8 & ~isempty(varargin{6})
tracename=varargin{6};
else
tracename='B';
end
[a,e]=rdaAUG_eff(shot,ppftype,tracename,shotfile_exp_eff,timerange);
trace=a;
trace.dim=[{[starti:endi]'} ; {trace.t}];
trace.x=trace.dim{1};
trace.dimunits=[{'channels'} ; {'time [s]'}];
trace.units='W/m^2';
trace.name=[num2str(shot) '/' ppftype '/' tracename];
% keep only nth points
trace.t=trace.t(1:nth:end);
trace.data=trace.data(:,1:nth:end);
trace.dim{2}=trace.t;
% calculating intersection of the view lines with magnetics axis
if strcmp(data_type_eff_noext,'sxR')
if nargin>=5 & ~isempty(varargin{3})
zmag=varargin{3};
else
zmag=loadAUGdata(shot,'zmag');
end
zmageff=interp1(zmag.t,zmag.data,trace.t);
if strcmp(tracename,'B')
[R_B, Z_B, ang_B,Rsxr]=sxrbgeometry(zmageff);
elseif strcmp(tracename,'A')
[R_A, Z_A, ang_A,Rsxr]=sxrageometry(zmageff);
else
disp(['sxr camera: ' tracename ' not set yet for calculating R projection'])
radius.data=Rsxr;
radius.t=trace.t;
varargout{1}={radius};
trace.R=radius;
end
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case {'sxb', 'sxf', 'ssx', 'ssx_g', 'ssx_h', 'ssx_i', 'ssx_j'}
% LOAD MULTI CHANNEL DATA SXB/J_0xx (or other than J camera if specified in varargin{8})
% load AUG soft x-ray data
if nargin>=3 & ~isempty(varargin{1})
% chords to be loaded
starti=varargin{1}(1);
endi=varargin{1}(2);
elseif strcmp(AUGkeywrdcase{index}(1:3),'ssx')
starti=1;
endi=60;
else
starti=52;
endi=54;
end
if nargin>=4 & ~isempty(varargin{2})
% chords already loaded, 1=not loaded (to load), 0=loaded
status=varargin{2};
else
status=ones(endi-starti+1,1);
end
if nargin>=6 & ~isempty(varargin{4})
timerange=varargin{4};
else
timerange=[0 10];
end
if nargin>=7 & ~isempty(varargin{5})
nth=varargin{5};
else
nth=13;
end
if nargin>=8 & ~isempty(varargin{6})
tracename=varargin{6};
else
if strcmp(AUGkeywrdcase{index},'sxb')
tracename='J';
elseif strcmp(AUGkeywrdcase{index},'sxf')
tracename='I';
elseif strcmp(AUGkeywrdcase{index},'ssx')
tracename='G';
elseif strcmp(AUGkeywrdcase{index}(1:4),'ssx_')
tracename=upper(AUGkeywrdcase{index}(5));
else
disp('should not be here, ask O. Sauter');
end
end
trace.t=[];
trace.x=[];
iok=0;
for ichord=starti:endi
tracename_eff = [tracename '_' num2str(ichord,'%.3d')];
try
[a,e]=rdaAUG_eff(shot,ppftype,tracename_eff,shotfile_exp_eff,timerange);
catch
a = [];
end
if isempty(a) || e~=0
if ~exist('trace_all')
trace_all = struct([]);
else
end
else
if iok==0
trace_all = a;
trace_all = rmfield(trace_all,[{'value'},{'data'}]);
iok = iok+1;
else
iok = iok+1;
end
trace_all.value(ichord,:) = a.value;
trace_all.data(ichord,:) = a.data;
end
end
if ~isempty(trace_all)
trace_all.dim=[{[starti:endi]'} ; {trace.t}];
trace = trace_all;
trace.x=trace.dim{1};
trace.dimunits=[{'channels'} ; {'time [s]'}];
trace.units='W/m^2';
trace.name=[num2str(shot) '/' ppftype '/' tracename];
% keep only nth points
trace.t=trace.t(1:nth:end);
trace.data=trace.data(:,1:nth:end);
trace.dim{2}=trace.t;
trace.value=trace.value(:,1:nth:end);
trace.time_aug.value=trace.time_aug.value(1:nth:end);
else
trace.data = [];
trace.dim = [];
trace.dimunits = [];
trace.x = [];
trace.t = [];
trace.units = [];
trace.name=[num2str(shot) '/' ppftype '/' tracename];
end
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
case {'ece','eced','ece_rho','eced_rho'}
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
% LOAD MULTI CHANNEL DATA
% load AUG ece data
shotfile_exp_eff = AUGexplocation{index};
if nargin>=3 & ~isempty(varargin{1})
starti=varargin{1}(1);
endi=varargin{1}(2);
else
starti=1;
endi=30;
end
if nargin>=4 & ~isempty(varargin{2})
status=varargin{2};
else
status=ones(endi-starti+1,1);
end
if nargin>=6 & ~isempty(varargin{4})
timerange=varargin{4};
else
timerange=[0 10];
end
if nargin>=7 & ~isempty(varargin{5})
nth=varargin{5};
else
nth=1;
end
trace.t=[];
trace.x=[];
ppftype=AUGsiglocation{1,index};
tracename=AUGsiglocation{2,index};
[a,e]=rdaAUG_eff(shot,ppftype,tracename,shotfile_exp_eff,timerange);
starti=1;
endi=size(a.value,1);
trace=a;
trace.dim=[{[starti:endi]'} ; {trace.t}];
trace.x=trace.dim{1};
trace.dimunits=[{'channels'} ; {'time [s]'}];
trace.units='W/m^2';
trace.name=[num2str(shot) '/' ppftype '/' tracename];
% keep only nth points
trace.t=trace.t(1:nth:end);
trace.data=trace.data(:,1:nth:end);
trace.dim{2}=trace.t;
% get R
[aR,e]=rdaAUG_eff(shot,ppftype,'R-A',shotfile_exp_eff,timerange);
[aZ,e]=rdaAUG_eff(shot,ppftype,'z-A',shotfile_exp_eff,timerange);
domatchRtime=0;
if domatchRtime
% interpolate R structure on ece data time array, to ease plot vs R
for i=starti:endi
radius.data(i,:) = interp1(aR.t,aR.data(i,:),trace.t);
zheight.data(i,:) = interp1(aZ.t,aZ.data(i,:),trace.t);
end
radius.t=trace.t;
zheight.t=trace.t;
else
radius.data = aR.data;
radius.t=aR.t;
zheight.data = aZ.data;
zheight.t=aR.t;
end
ij=find(trace.data==0);
trace.data(ij)=NaN;
varargout{1}={radius};
trace.R=radius;
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
trace.Z=zheight;
if strcmp(AUGkeywrdcase{index},'ece_rho') || strcmp(AUGkeywrdcase{index},'eced_rho')
equil=gdat(shot,'equil',0);
inb_chord_ece=size(trace.R.data,1);
inb_time_ece=size(trace.R.data,2);
psi_out = NaN*ones(inb_chord_ece,inb_time_ece);
rhopsinorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
rhotornorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
rhovolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-1)];
iok=find(aR.data(:,1)>0);
for itequil=1:length(time_equil)-1
rr=equil.Rmesh(:,itequil);
zz=equil.Zmesh(:,itequil);
psirz_in = equil.psi2D(:,:,itequil);
it_ece_inequil = find(trace.R.t>=time_equil(itequil) & trace.R.t<=time_equil(itequil+1));
if ~isempty(it_ece_inequil)
rout=trace.R.data(iok,it_ece_inequil);
ijok=find(~isnan(rout));
if ~isempty(ijok)
zout=trace.Z.data(iok,it_ece_inequil);
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
psi_out(iok,it_ece_inequil) = reshape(psi_at_routzout,length(iok),length(it_ece_inequil));
rhopsinorm_out(iok,it_ece_inequil) = sqrt((psi_out(iok,it_ece_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
for it_cx=1:length(it_ece_inequil)
rhotornorm_out(iok,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out(iok,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
rhovolnorm_out(iok,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out(iok,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
end
end
trace.rhos.psi_on_rztime = psi_out;
trace.rhos.rhopsinorm_on_rztime = rhopsinorm_out;
trace.rhos.rhotornorm_on_rztime = rhotornorm_out;
trace.rhos.rhovolnorm_on_rztime = rhovolnorm_out;
trace.rhos.t = trace.R.t;
end
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
% LOAD CEZ data
%
shotfile_exp_eff = AUGexplocation{index};
trace.t=[];
trace.x=[];
ppftype='CEZ';
[a,e]=rdaAUG_eff(shot,ppftype,'vrot',shotfile_exp_eff);
if isempty(a) || e~=0
cxrs = struct([]);
else
cxrs.vrot = a;
rmfield(cxrs.vrot,'x'); rmfield(cxrs.vrot,'t'); rmfield(cxrs.vrot,'data');
cxrs.x = a.x;
cxrs.data = a.data;
cxrs.t = a.t;
cxrs.x = a.x;
end
if ~isempty(cxrs)
cxrs.dim=[{cxrs.x} ; {cxrs.t}];
cxrs.dimunits=[{'chord'} ; {'time [s]'}];
cxrs.name=[num2str(shot) '/' ppftype '/vrot;Ti;Ti_c'];
[aerr,e]=rdaAUG_eff(shot,ppftype,'err_vrot',shotfile_exp_eff);
cxrs.vrot.error = aerr.value;
[r_time]=sf2ab(ppftype,shot,'R_time','-exp',shotfile_exp_eff);
cxrs.r_time = r_time.value{1};
[z_time]=sf2ab(ppftype,shot,'z_time','-exp',shotfile_exp_eff);
cxrs.z_time = z_time.value{1};
[a,e]=rdaAUG_eff(shot,ppftype,'Ti',shotfile_exp_eff);
[aerr,e]=rdaAUG_eff(shot,ppftype,'err_Ti',shotfile_exp_eff);
cxrs.ti = a;
cxrs.ti.error = aerr.value;
[a,e]=rdaAUG_eff(shot,ppftype,'Ti_c',shotfile_exp_eff);
[aerr,e]=rdaAUG_eff(shot,ppftype,'err_Ti_c',shotfile_exp_eff);
cxrs.ti_c = a;
cxrs.ti_c.error = aerr.value;
%
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
equil=gdat(shot,'equil',0);
inb_chord_cxrs=size(cxrs.r_time,1);
inb_time_cxrs=size(cxrs.r_time,2);
psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhopsinorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-1)];
iok=find(cxrs.r_time(:,1)>0);
for itequil=1:length(time_equil)-1
rr=equil.Rmesh(:,itequil);
zz=equil.Zmesh(:,itequil);
psirz_in = equil.psi2D(:,:,itequil);
it_cxrs_inequil = find(cxrs.t>=time_equil(itequil) & cxrs.t<=time_equil(itequil+1));
if ~isempty(it_cxrs_inequil)
rout=cxrs.r_time(iok,it_cxrs_inequil);
zout=cxrs.z_time(iok,it_cxrs_inequil);
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil));
rhopsinorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
for it_cx=1:length(it_cxrs_inequil)
rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
end
cxrs.psi_on_rztime = psi_out;
cxrs.rhopsinorm_on_rztime = rhopsinorm_out;
cxrs.rhotornorm_on_rztime = rhotornorm_out;
cxrs.rhovolnorm_on_rztime = rhovolnorm_out;
end
trace = cxrs;
else
trace.data = [];
trace.dim = [];
trace.dimunits = [];
trace.x = [];
trace.t = [];
trace.units = [];
trace.name=[num2str(shot) '/' ppftype '/' 'vrot, Ti_c,...'];
end
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
case {'neterho'} % 'nerho', 'terho',
% LOAD VTA data
% Vertical Thomson core and edge
shotfile_exp_eff = AUGexplocation{index};
trace.t=[];
trace.x=[];
ppftype='VTA';
% if strcmp(AUGkeywrdcase{index},'terho')
[a,e]=rdaAUG_eff(shot,ppftype,'Te_c',shotfile_exp_eff);
if isempty(a) || e~=0
thomson = struct([]);
else
thomson.te_core = a;
rmfield(thomson.te_core,'x'); rmfield(thomson.te_core,'t'); rmfield(thomson.te_core,'data');
thomson.x = a.x;
thomson.data = a.data;
thomson.t = a.t;
thomson.time_core = a.t;
thomson.x = a.x;
end
if ~isempty(thomson)
thomson.dim=[{thomson.x} ; {thomson.t}];
thomson.dimunits=[{'points'} ; {'time [s]'}];
thomson.name=[num2str(shot) '/' ppftype '/te_c;te_e;ne_c;ne_e'];
[alow,e]=rdaAUG_eff(shot,ppftype,'Telow_c',shotfile_exp_eff);
[aup,e]=rdaAUG_eff(shot,ppftype,'Teupp_c',shotfile_exp_eff);
thomson.te_core.error = max(aup.value-thomson.te_core.value,thomson.te_core.value-alow.value);
%
[a,e]=rdaAUG_eff(shot,ppftype,'Ne_c',shotfile_exp_eff);
thomson.ne_core = a;
rmfield(thomson.ne_core,'x'); rmfield(thomson.ne_core,'t'); rmfield(thomson.ne_core,'data');
[alow,e]=rdaAUG_eff(shot,ppftype,'Nelow_c',shotfile_exp_eff);
[aup,e]=rdaAUG_eff(shot,ppftype,'Neupp_c',shotfile_exp_eff);
thomson.ne_core.error = max(aup.value-thomson.ne_core.value,thomson.ne_core.value-alow.value);
%
[r_time]=rdaAUG_eff(shot,ppftype,'R_core',shotfile_exp_eff);
thomson.r_core_time = ones(size(thomson.te_core.value,1),1)*r_time.value;
[z_time]=rdaAUG_eff(shot,ppftype,'Z_core',shotfile_exp_eff);
thomson.z_core_time = z_time.value' * ones(1,size(thomson.te_core.value,2));
%
[a,e]=rdaAUG_eff(shot,ppftype,'Te_e',shotfile_exp_eff);
thomson.te_edge = a;
thomson.time_edge = a.t;
rmfield(thomson.te_edge,'x'); rmfield(thomson.te_edge,'t'); rmfield(thomson.te_edge,'data');
[alow,e]=rdaAUG_eff(shot,ppftype,'Telow_e',shotfile_exp_eff);
[aup,e]=rdaAUG_eff(shot,ppftype,'Teupp_e',shotfile_exp_eff);
thomson.te_edge.error = max(aup.value-thomson.te_edge.value,thomson.te_edge.value-alow.value);
%
[a,e]=rdaAUG_eff(shot,ppftype,'Ne_e',shotfile_exp_eff);
thomson.ne_edge = a;
rmfield(thomson.ne_edge,'x'); rmfield(thomson.ne_edge,'t'); rmfield(thomson.ne_edge,'data');
[alow,e]=rdaAUG_eff(shot,ppftype,'Nelow_e',shotfile_exp_eff);
[aup,e]=rdaAUG_eff(shot,ppftype,'Neupp_e',shotfile_exp_eff);
thomson.ne_edge.error = max(aup.value-thomson.ne_edge.value,thomson.ne_edge.value-alow.value);
%
[r_time]=rdaAUG_eff(shot,ppftype,'R_edge',shotfile_exp_eff);
thomson.r_edge_time = ones(size(thomson.te_edge.value,1),1)*r_time.value;
[z_time]=rdaAUG_eff(shot,ppftype,'Z_edge',shotfile_exp_eff);
thomson.z_edge_time = z_time.value' * ones(1,size(thomson.te_edge.value,2));
%
if strcmp(AUGkeywrdcase{index},'neterho')
equil=gdat(shot,'equil',0);
% core
inb_chord_thomson_core=size(thomson.r_core_time,1);
inb_time_thomson_core=size(thomson.r_core_time,2);
psi_out_core = NaN*ones(inb_chord_thomson_core,inb_time_thomson_core);
rhopsinorm_out_core = NaN*ones(inb_chord_thomson_core,inb_time_thomson_core);
rhotornorm_out_core = NaN*ones(inb_chord_thomson_core,inb_time_thomson_core);
rhovolnorm_out_core = NaN*ones(inb_chord_thomson_core,inb_time_thomson_core);
% edge
inb_chord_thomson_edge=size(thomson.r_edge_time,1);
inb_time_thomson_edge=size(thomson.r_edge_time,2);
psi_out_edge = NaN*ones(inb_chord_thomson_edge,inb_time_thomson_edge);
rhopsinorm_out_edge = NaN*ones(inb_chord_thomson_edge,inb_time_thomson_edge);
rhotornorm_out_edge = NaN*ones(inb_chord_thomson_edge,inb_time_thomson_edge);
rhovolnorm_out_edge = NaN*ones(inb_chord_thomson_edge,inb_time_thomson_edge);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-1)];
for itequil=1:length(time_equil)-1
rr=equil.Rmesh(:,itequil);
zz=equil.Zmesh(:,itequil);
psirz_in = equil.psi2D(:,:,itequil);
it_thomson_core_inequil = find(thomson.time_core>=time_equil(itequil) & thomson.time_core<=time_equil(itequil+1));
if ~isempty(it_thomson_core_inequil)
rout_core=thomson.r_core_time(:,it_thomson_core_inequil);
zout_core=thomson.z_core_time(:,it_thomson_core_inequil);
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core);
psi_out_core(:,it_thomson_core_inequil) = reshape(psi_at_routzout,inb_chord_thomson_core,length(it_thomson_core_inequil));
rhopsinorm_out_core(:,it_thomson_core_inequil) = sqrt((psi_out_core(:,it_thomson_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
for it_cx=1:length(it_thomson_core_inequil)
rhotornorm_out_core(:,it_thomson_core_inequil(it_cx)) = ...
interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out_core(:,it_thomson_core_inequil(it_cx)),-3,[2 2],[0 1]);
rhovolnorm_out_core(:,it_thomson_core_inequil(it_cx)) = ...
interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out_core(:,it_thomson_core_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
% edge
it_thomson_edge_inequil = find(thomson.time_edge>=time_equil(itequil) & thomson.time_edge<=time_equil(itequil+1));
if ~isempty(it_thomson_edge_inequil)
rout_edge=thomson.r_edge_time(:,it_thomson_edge_inequil);
zout_edge=thomson.z_edge_time(:,it_thomson_edge_inequil);
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_edge,zout_edge);
psi_out_edge(:,it_thomson_edge_inequil) = reshape(psi_at_routzout,inb_chord_thomson_edge,length(it_thomson_edge_inequil));
rhopsinorm_out_edge(:,it_thomson_edge_inequil) = sqrt((psi_out_edge(:,it_thomson_edge_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
for it_cx=1:length(it_thomson_edge_inequil)
rhotornorm_out_edge(:,it_thomson_edge_inequil(it_cx)) = ...
interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out_edge(:,it_thomson_edge_inequil(it_cx)),-3,[2 2],[0 1]);
rhovolnorm_out_edge(:,it_thomson_edge_inequil(it_cx)) = ...
interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out_edge(:,it_thomson_edge_inequil(it_cx)),-3,[2 2],[0 1]);