add rhovol IOH, etc

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@1895 d63d8f72-b253-0410-a779-e742ad2e26cf

TCV/loadTCVdata.m

  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 boundary vs time
-% 'zcont' =  Z of plama boundary vs time
-% 'vol' =  volume of flux surfaces
-% 'qrho' =  q profile on rho mesh
-% 'q95' =  q95 vs time
-% 'kappa', 'elon' =  edge elongation vs time
-% 'delta', 'triang' =  edge averaged triangularity vs time
-% 'deltatop', 'triangtop' =  edge upper (top) triangularity vs time
-% 'deltabot', 'triangbot' =  edge lower (bottom) triangularity vs time
+% list of data_type currently available (when [_2,_3] is added, means can add _i to get Liuqe i):
+% 'Ip'[_2,_3]   =  current
+% 'zmag'[_2,_3] =  vertical position of the center of the plasma (magnetic axis)
+% 'rmag'[_2,_3] =  radial position of the center of the plasma
+% 'rcont'[_2,_3] =  R of plama boundary vs time
+% 'zcont'[_2,_3] =  Z of plama boundary vs time
+% 'vol'[_2,_3] =  volume of flux surfaces
+% 'rhovol'[_2,_3] =  sqrt(V(:,t)/V(edge,t)), normalised rho variable based on volume of flux surfaces
+% 'qrho'[_2,_3] =  q profile on rho mesh
+% 'q95'[_2,_3] =  q95 vs time
+% 'kappa', 'elon'[_2,_3] =  edge elongation vs time
+% 'delta', 'triang'[_2,_3] =  edge averaged triangularity vs time
+% 'deltatop', 'triangtop'[_2,_3] =  edge upper (top) triangularity vs time
+% 'deltabot', 'triangbot'[_2,_3] =  edge lower (bottom) triangularity vs time
 % '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
+% '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)
 % 'neft' =  ne fitted from data on rho mesh (from proffit.local_time:neft_abs)
@@ -28,6 +29,8 @@
 % 'sxr'  =  soft x-ray emission
 % 'sxR'  =  soft x-ray emission with varargout{1} option (requires varargin{4}!)
 % 'MPX'  =  soft x-ray from wire chambers
+% 'IOH'  =  current in ohmic transformer coils IOH_1
+% 'vloop'  =  loop voltage
 %
 % 'xx_2 or xx_3' for Liuqe2 or 3: same as above for xx related to equilibrium
 %
@@ -38,11 +41,12 @@
 % Definition of varargin depends on data_type:
 %
 % data_type=sxr or ece:
-%                  varargin{1}:  channel status 1= unread yet, 0= read
+%                  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{2}:  [i1 i2] : if not empty, assumes need many chords from i1 to i2
 %                  varargin{3}:  zmag for varargout{1} computation
+%                       (can have more inputs for AUG, but not used here)
 %
 % OUTPUT:
 %
@@ -107,6 +111,9 @@ end
 if ~isempty(strmatch(data_type_eff_noext,[{'VOL'} {'volume'}],'exact'))
   data_type_eff_noext='vol';
 end
+if ~isempty(strmatch(data_type_eff_noext,[{'RHOVOL'}],'exact'))
+  data_type_eff_noext='rhovol';
+end
 if ~isempty(strmatch(data_type_eff_noext,[{'q_95'} {'Q95'}],'exact'))
   data_type_eff_noext='q95';
 end
@@ -128,6 +135,15 @@ end
 if ~isempty(strmatch(data_type_eff_noext,[{'Zmag'}],'exact'))
   data_type_eff_noext='zmag';
 end
+if ~isempty(strmatch(data_type_eff_noext,[{'MPX'}],'exact'))
+  data_type_eff_noext='MPX';
+end
+if ~isempty(strmatch(data_type_eff_noext,[{'ioh'} {'Ioh'} {'iot'} {'IOT'}],'exact'))
+  data_type_eff_noext='IOH';
+end
+if ~isempty(strmatch(data_type_eff_noext,[{'Vloop'} {'Vsurf'} {'vsurf'}],'exact'))
+  data_type_eff_noext='vloop';
+end
 
 % some defaults
 
@@ -135,20 +151,21 @@ end
 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'}];
+      {'\results::thomson:psiscatvol'} {'\results::thomson:psi_max'}];
 
 % all keywords and corresponding case to run below
-TCVkeywrdall=[{'Ip'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'qrho'} {'q95'} {'kappa'} ...
+TCVkeywrdall=[{'Ip'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'rhovol'} {'qrho'} {'q95'} {'kappa'} ...
       {'delta'} {'deltatop'} {'deltabot'} {'neint'} ...
       {'ne'} {'te'} {'nerho'} {'terho'} {'profnerho'} {'profterho'} ...
       {'neft'} {'teft'} {'neftav'} {'teftav'}  ...
-      {'sxr'} {'sxR'} {'ece'}];
+      {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'}];
 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.irhovol=strmatch('rhovol',TCVkeywrdall,'exact');
 TCVsig.iqrho=strmatch('qrho',TCVkeywrdall,'exact');
 TCVsig.iq95=strmatch('q95',TCVkeywrdall,'exact');
 TCVsig.ikappa=strmatch('kappa',TCVkeywrdall,'exact');
@@ -169,12 +186,17 @@ TCVsig.iteftav=strmatch('teftav',TCVkeywrdall,'exact');
 TCVsig.isxr=strmatch('sxr',TCVkeywrdall,'exact');
 TCVsig.isxR=strmatch('sxR',TCVkeywrdall,'exact');
 TCVsig.iece=strmatch('ece',TCVkeywrdall,'exact');
+TCVsig.iMPX=strmatch('MPX',TCVkeywrdall,'exact');
+TCVsig.iIOH=strmatch('IOH',TCVkeywrdall,'exact');
+TCVsig.ivloop=strmatch('vloop',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.ivol)=TCVkeywrdall(TCVsig.ivol); % special as nodes _2 or _3 not existing with psitbx
+TCVkeywrdcase(TCVsig.irhovol)=TCVkeywrdall(TCVsig.irhovol); % idem vol
 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
@@ -182,6 +204,9 @@ TCVkeywrdcase(TCVsig.iterho)=TCVkeywrdall(TCVsig.iterho); % idem
 TCVkeywrdcase(TCVsig.isxr)=TCVkeywrdall(TCVsig.isxr);
 TCVkeywrdcase(TCVsig.isxR)=TCVkeywrdall(TCVsig.isxR);
 TCVkeywrdcase(TCVsig.iece)=TCVkeywrdall(TCVsig.iece);
+TCVkeywrdcase(TCVsig.iMPX)=TCVkeywrdall(TCVsig.iMPX);
+TCVkeywrdcase(TCVsig.iIOH)=TCVkeywrdall(TCVsig.iIOH);
+TCVkeywrdcase(TCVsig.ivloop)=TCVkeywrdall(TCVsig.ivloop);
 
 % 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
@@ -196,7 +221,6 @@ 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.iq95)={'\results::q_95'};
 TCVsiglocation(TCVsig.ikappa)={'\results::kappa_edge'};
 TCVsiglocation(TCVsig.idelta)={'\results::delta_edge'};
@@ -225,7 +249,7 @@ if strcmp(data_type_eff(1:1),'\')
     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'])
+%    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;
@@ -253,32 +277,16 @@ else
 end
 disp(['loading' ' ' data_type_eff_noext ' from TCV shot #' num2str(shot)]);
 disp(['case ' TCVkeywrdcase{index}])
-if i_23==2 & isempty(strmatch(TCVsiglocation(index),liuqe23,'exact'))
+if i_23==2 & isempty(strmatch(TCVsiglocation(index),liuqe23,'exact')) & strcmp(TCVkeywrdcase{index},'simpletdi')
   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'])
+%  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}
 
@@ -290,8 +298,10 @@ switch TCVkeywrdcase{index}
     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)])
+    ij=findstr(nodenameeff,'[');
+    if isempty(ij); ij=length(nodenameeff)+1; end
+    if eval(['~mdsdata(''node_exists("\' nodenameeff(1:ij-1) '")'')'])
+      disp(['node ' nodenameeff(1:ij-1) ' does not exist for shot = ' num2str(shot)])
       return
 % $$$     elseif eval(['mdsdata(''getnci("\' nodenameeff ':foo","length")'')==0'])
 % $$$       disp(['no data for node ' nodenameeff ' for shot = ' num2str(shot)])
@@ -304,15 +314,21 @@ switch TCVkeywrdcase{index}
       return
     end
     trace.data=tracetdi.data;
+
     if TCVsigtimeindx(index)>0 | length(tracetdi.dim)==1
       trace.t=tracetdi.dim{max(1,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
+      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
     if length(tracetdi.dim)==2
@@ -369,7 +385,7 @@ switch TCVkeywrdcase{index}
     trace.dimunits=[{'Z [m]'} ; {'time [s]'}];
     trace.x=z;
     trace.t=time;
-    mdsclose
+    mdsclose('mdsip')
 
   %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
   case {'nerho','terho'}
@@ -390,6 +406,7 @@ 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})
@@ -415,9 +432,70 @@ switch TCVkeywrdcase{index}
     trace.name=[num2str(shot) ';' nodenameeff];
     mdsclose
 
+  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
+  case {'vol'}
+    % vol from psitbx
+    mdsopen(shot);
+    nodenameeff=['\results::psitbx:vol'];
+    tracetdi=tdi(nodenameeff);
+    if i_23==2
+      ['LIUQE' liuqe_ext(2:2)]
+      psi=psitbxtcv(shot,'+0',['LIUQE' liuqe_ext(2:2)]);
+      fsd=psitbxp2p(psi,'FS');
+      fsg=psitbxfsg(fsd);
+      tracetdi.data = fsg.vol.x;
+      tracetdi.dim{1}=fsg.vol.grid.x{:}';
+      tracetdi.dim{2}=fsg.vol.grid.t';
+    end
+    trace.data=tracetdi.data;
+    trace.x=tracetdi.dim{1};
+    trace.t=tracetdi.dim{2};
+    trace.dim=[{trace.x};{trace.t}];
+    trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
+    trace.name=[num2str(shot) ';' nodenameeff liuqe_ext];
+    mdsclose
+
+  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
+  case {'rhovol'}
+    % vol from psitbx
+    mdsopen(shot);
+    nodenameeff=['\results::psitbx:vol'];
+    tracetdi=tdi(nodenameeff);
+    if i_23==2
+      ['LIUQE' liuqe_ext(2:2)]
+      psi=psitbxtcv(shot,'+0',['LIUQE' liuqe_ext(2:2)]);
+      fsd=psitbxp2p(psi,'FS');
+      fsg=psitbxfsg(fsd);
+      tracetdi.data = fsg.vol.x;
+      tracetdi.dim{1}=fsg.vol.grid.x{:}';
+      tracetdi.dim{2}=fsg.vol.grid.t';
+    end
+    trace.data=tracetdi.data;
+    for i=1:size(tracetdi.data,2)
+      trace.data(:,i)=sqrt(tracetdi.data(:,i)./tracetdi.data(end,i));
+    end
+    trace.x=tracetdi.dim{1};
+    trace.t=tracetdi.dim{2};
+    trace.dim=[{trace.x};{trace.t}];
+    trace.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
+    trace.name=[num2str(shot) '; sqrt(V/Va) from ' nodenameeff liuqe_ext];
+    mdsclose
+
   %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
   case {'sxr','sxR'}
     %  load TCV soft x-ray data
+    if  nargin>=3 & ~isempty(varargin{1})
+      i1=varargin{1}(1);
+      i2=varargin{1}(2);
+    else
+      i1=1;
+      i2=20;
+    end
+    if  nargin>=4 & ~isempty(varargin{2})
+      status=varargin{2};
+    else
+      status=ones(i2-i1+1,1);
+    end
 
     % 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
@@ -425,7 +503,9 @@ switch TCVkeywrdcase{index}
       [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)
+        if nargin>=5 & ~isempty(varargin{3})
+          zmag=varargin{3};
+        else
           zmag=loadTCVdata(shot,['zmag' liuqe_ext]);
         end
         t_1=zmag.t(1);
@@ -495,6 +575,7 @@ switch TCVkeywrdcase{index}
     t_2=zmag.t(end);
     if ~isempty(find(status == 1))
       mdsopen(shot);
+      keyboard
       signal=get_mds_mio('MPX',[t_1 t_2]);
       mdsclose
       trace.data=signal.data;
@@ -509,6 +590,36 @@ switch TCVkeywrdcase{index}
     varargout{1}={VsxrTCVradius(zmag.data,xchord,ychord)};
     error=0;
 
+  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
+  case 'IOH'
+    %  Ohmic transformer current
+    mdsopen(shot);
+    nodenameeff=[{'\magnetics::ipol[*,$1]'} {'OH_001'}];
+    tracetdi=tdi(nodenameeff{:});
+    trace.data=tracetdi.data;
+    trace.x=[];
+    trace.t=tracetdi.dim{1};
+    trace.dim=tracetdi.dim;
+    trace.dimunits={'time [s]'};
+    trace.name=[num2str(shot) ';' nodenameeff{1} ',' nodenameeff{2}];
+    mdsclose
+    error=0;
+    
+  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
+  case 'vloop'
+    %  Loop voltage
+    mdsopen(shot);
+    nodenameeff=[{'\magnetics::vloop[*,$1]'} {'001'}];
+    tracetdi=tdi(nodenameeff{:});
+    trace.data=tracetdi.data;
+    trace.x=[];
+    trace.t=tracetdi.dim{1};
+    trace.dim=tracetdi.dim;
+    trace.dimunits={'time [s]'};
+    trace.name=[num2str(shot) ';' nodenameeff{1} ',' nodenameeff{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']);