Skip to content
Snippets Groups Projects
Commit 56ff35d7 authored by Olivier Sauter's avatar Olivier Sauter
Browse files

add AUG

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@1865 d63d8f72-b253-0410-a779-e742ad2e26cf
parent 7f3cb9f3
No related branches found
No related tags found
No related merge requests found
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}!)
%
% gdat(15133,'MAG/Ipi',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'
% data_type='POT/ELMa-Han', 'MOD/OddNAmp', 'MOD/EvenNAmp', 'TOT/PNBI_TOT'
%
% Meaning 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
%
% 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');
%
varargout={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;
if size(data_type_eff,1)==1
i=findstr('/',data_type_eff);
if length(i)>=1
% assumes given a la 'MAG/Ipi'
data_type_eff=[{data_type_eff(1:i(1)-1)} ; {data_type_eff(i(1)+1:end)}];
end
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'} {'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,[{'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,[{'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,[{'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,[{'Rcont'}],'exact'))
data_type_eff_noext='rcont';
end
if ~isempty(strmatch(data_type_eff_noext,[{'Zcont'}],'exact'))
data_type_eff_noext='zcont';
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'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'qrho'} {'q95'} {'kappa'} ...
{'delta'} {'deltatop'} {'deltabot'} {'neint'} ...
{'ne'} {'te'} {'nerho'} {'terho'} ...
{'sxr'} {'sxR'} {'ece'}];
AUGsig.iip=strmatch('Ip',AUGkeywrdall,'exact');
AUGsig.izmag=strmatch('zmag',AUGkeywrdall,'exact');
AUGsig.irmag=strmatch('rmag',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.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.isxr=strmatch('sxr',AUGkeywrdall,'exact');
AUGsig.isxR=strmatch('sxR',AUGkeywrdall,'exact');
AUGsig.iece=strmatch('ece',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.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.iterho)=AUGkeywrdall(AUGsig.iterho); % idem
AUGkeywrdcase(AUGsig.isxr)=AUGkeywrdall(AUGsig.isxr);
AUGkeywrdcase(AUGsig.isxR)=AUGkeywrdall(AUGsig.isxR);
%AUGkeywrdcase(AUGsig.iece)=AUGkeywrdall(AUGsig.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
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'; 'Ipi'};
AUGsiglocation(:,AUGsig.izmag)={''; ''};
AUGsiglocation(:,AUGsig.irmag)={''; ''};
AUGsiglocation(:,AUGsig.ircont)={'' ; ''}; AUGsigtimeindx(AUGsig.ircont)=2;
AUGsiglocation(:,AUGsig.izcont)={'' ; ''}; AUGsigtimeindx(AUGsig.izcont)=2;
AUGsiglocation(:,AUGsig.ivol)={''; ''};
AUGsiglocation(:,AUGsig.iq95)={''; ''};
AUGsiglocation(:,AUGsig.ikappa)={''; ''};
AUGsiglocation(:,AUGsig.ideltatop)={''; ''};
AUGsiglocation(:,AUGsig.ideltabot)={''; ''};
AUGsiglocation(:,AUGsig.ineint)={'DCN'; 'H-1'};
% 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}];
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~='''');
tracename=tracename(ij);
[a,e]=rdaAUG_eff(shot,ppftype,tracename);
% switch tracename
% special cases if traces do not exist for some shot or other
% end
trace=a;
clear error
error=e;
if 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?')
trace.dim=[{trace.x} ; {trace.t} ; {d}];
trace.dimunits=[{'R [m] or rho=sqrt(psi_norm)'} ; {'time [s]'} ; {'d'}];
trace.d=d;
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
trace.t=[];
trace.x=[];
ppftype='SXR';
tracename='B';
[a,e]=rdaAUG_eff(shot,ppftype,tracename);
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];
% 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);
% for i=starti:endi
% radius.data(i,:)=2.848 + (2.172-zmageff') .* tan(-4.5/180.*3.14159 - atan2(0.99.*(i-18),35.31));
% end
% radius.t=t;
varargout{1}={struct([])};
trace.R=[];
% end
otherwise
disp('case not yet defined')
end
function [adata,error]=rdaAUG_eff(shot,diagname,sigtype,varargin);
%
% gets data using SF2ML or mdsplus (when mdsplus will be available)
% 1D arrays: assumes dimension is time
% 2D arrays: assumes data vs (x,time)
% 3D arrays: assumes data vs (x,time,hsig) (for mdsplus)
%
% varargin{1}: time interval or timevalue, will get data closest to that time or within that time interval
%
% examples:
% [data,error]=rdaAUG_eff(15133,'MAG','Ipi');
% [data,error]=rdaAUG_eff(15133,'MAG','Ipi',[1 5]);
%
% set global variable: usemdsplus to decide if SF2ML or mdsplus is used:
% >> global usemdsplus
% >> usemdsplus=1 % means use mds to get data
% >> usemdsplus=0 % means use SF2ML (default if not defined)
% if ~exist('usemdsplus'); usemdsplus=0; end
%
global usemdsplus
if isempty(usemdsplus); usemdsplus=0; end
error=1;
time_int=[];
if nargin>=4 & ~isempty(varargin{1})
time_int=varargin{1};
end
if usemdsplus
else
% use SF2ML
if isempty(time_int)
adata=SF2ML(sigtype,diagname,shot);
else
adata=SF2ML(sigtype,diagname,shot,time_int(1),time_int(end));
end
if size(adata.y,1)==1
adata.x=[];
adata.y=adata.y';
else
adata.x=[1:size(adata.y,1)];
end
adata.data=adata.y;
adata.t=adata.t';
% % transpose data as output in C format, reversed from Fortran and matlab standard
% ss=size(a);
% nbofdim=length(ss);
% if ss(end)==1; nbofdim=nbofdim-1; end
% nbofdim=max(nbofdim,1);
% if nbofdim==1
% data=a;
% else
% data=a';
% end
end
error=0;
function [data,x,time,hsig,error]=rda_eff(shot,pftype,tracename,varargin);
%
% gets data using RDA or mdsplus
% 1D arrays: assumes dimension is time
% 2D arrays: assumes data vs (x,time)
% 3D arrays: assumes data vs (x,time,hsig) (for mdsplus)
%
% varargin{1}: time interval or timevalue, will get data closest to that time or within that time interval
% (DOES NOT WORK YET)
%
% examples:
% [data,x,time,hsig,error]=rda_eff(51994,'ppf','efit/xip');
% [data,x,time,hsig,error]=rda_eff(52206,'ppf','equi/rmji?uid=jetthg+seq=122');
%
% set global variable: usemdsplus to decide if RDA or mdsplus is used:
% >> global usemdsplus
% >> usemdsplus=1 % means use mds to get data (default if not defined)
% >> usemdsplus=0 % means use jetreaddata routine (RDA)
% if ~exist('usemdsplus'); usemdsplus=1; end
%
global usemdsplus
if isempty(usemdsplus); usemdsplus=1; end
time_int=[];
if nargin>=4 & ~isempty(varargin{1})
time_int=varargin{1};
end
if usemdsplus
% use mdsplus
if ~unix('test -d /home/duval/mdsplus')
addpath('/home/duval/mdsplus')
end
mdsconnect('mdsplus.jet.efda.org');
% defines trace to fetch
% after '?' specific details
separator='+';
mainseparator='?';
imaintrace=findstr(mainseparator,tracename);
if isempty(imaintrace)
maintrace=tracename;
uid=[];
seq=[];
diag=[];
else
maintrace=tracename(1:imaintrace-1);
rest=tracename(imaintrace+1:end);
% gets uid if any
iuid=findstr('uid=',rest);
if isempty(iuid)
uid=[];
else
ii=findstr(separator,rest(iuid:end));
if isempty(ii)
uid=rest(iuid+4:end);
else
uid=rest(iuid+4:iuid+ii(1)-2);
end
end
% gets seq if any
iseq=findstr('seq=',rest);
if isempty(iseq)
seq=[];
else
ii=findstr(separator,rest(iseq:end));
if isempty(ii)
seq=rest(iseq+4:end);
else
seq=rest(iseq+4:iseq+ii(1)-2);
end
end
% gets diag if any
idiag=findstr('diag=',rest);
if isempty(idiag)
diag=[];
else
ii=findstr(separator,rest(idiag:end));
if isempty(ii)
diag=rest(idiag+5:end);
else
diag=rest(idiag+5:idiag+ii(1)-2);
end
end
end
% fetch value
if ~isempty(uid)
eval(['u=mdsvalue(''_sig=ppfuid("' uid '")'');'])
end
traceeff=[pftype '/' maintrace];
if ~isempty(seq)
traceeff=[traceeff '/' num2str(seq)];
end
user=getenv('USER');
eval(['[data,error]=mdsvalue(''_rdaeff' user '=jet("' traceeff '",' num2str(shot) ')'');'])
hsig=[];
ss=size(data);
nbofdim=length(ss);
if ss(end)==1; nbofdim=nbofdim-1; end
nbofdim=max(nbofdim,1);
switch nbofdim
case 1
eval(['time=mdsvalue(''dim_of(_rdaeff' user ',0)'');']);
x=[];
case 2
eval(['x=mdsvalue(''dim_of(_rdaeff' user ',0)'');']);
eval(['time=mdsvalue(''dim_of(_rdaeff' user ',1)'');']);
case 3
eval(['x=mdsvalue(''dim_of(_rdaeff' user ',0)'');']);
eval(['time=mdsvalue(''dim_of(_rdaeff' user ',1)'');']);
disp('3rd dimension in hsig!!!!!!!!!!!!!!!!!!!!!!!!!')
eval(['hsig=mdsvalue(''dim_of(_rdaeff' user ',2)'');']);
otherwise
disp([' more than 3 dimensions for ' num2str(shot) ' ; ' pftype '/' tracename])
error('in rda_eff')
end
mdsdisconnect;
if ~unix('test -d /home/duval/mdsplus')
rmpath('/home/duval/mdsplus')
end
else
% use RDA
[a,time,x,hsig,error]=jetreaddata(['http://data.jet.uk/' pftype '/' num2str(shot) '/' tracename]);
% transpose data as output in C format, reversed from Fortran and matlab standard
ss=size(a);
nbofdim=length(ss);
if ss(end)==1; nbofdim=nbofdim-1; end
nbofdim=max(nbofdim,1);
if nbofdim==1
data=a;
else
data=a';
end
end
% to prevent problems when trace empty and time become string
if ischar(time)
time=[];
end
if ischar(x)
x=[];
end
if isempty(x) & ~isempty(data) & data==0
data=[];
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment