diff --git a/matlab/JET/jet_requests_mapping.m b/matlab/JET/jet_requests_mapping.m
index 4ed7576f49814d24d16e7af2c3f8655fe71c2749..664bf4be1a85890b5eb2aebb0e9b888cb37ad249 100644
--- a/matlab/JET/jet_requests_mapping.m
+++ b/matlab/JET/jet_requests_mapping.m
@@ -340,10 +340,20 @@ switch lower(data_request)
   mapping.timedim = 1;
   mapping.method = 'signal';
   mapping.expression = [{'PPF'},{'EFIT'},{'RMAG'}];
+ case {'rnt', 'neutron_rate'}
+  mapping.label = 'neutron\_rate';
+  mapping.timedim = 1;
+  mapping.method = 'signal';
+  mapping.expression = [{'PPF'},{'TIN'},{'RNT'}];
  case {'sxr', 'sxr_s40', 'sxr_htv'}
   mapping.timedim = 1;
   mapping.gdat_timedim = 2;
   mapping.method = 'switchcase';
+ case {'tbeo'}
+  mapping.label = 'TBEO';
+  mapping.timedim = 1;
+  mapping.method = 'signal';
+  mapping.expression = [{'PPF'},{'EDG8'},{'TBEO'}];
  case 'te'
   mapping.timedim = 2;
   mapping.label = 'Te';
diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index b4e524f427ecf5ff4a435439a6a1a24af9f9d9a5..ba6cab40f7388c591dccadd87cda1db222ad8fa4 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -312,10 +312,10 @@ end
 
 mapping_for_tcv = tcv_requests_mapping(data_request_eff,shot);
 gdat_data.label = mapping_for_tcv.label;
-% special treatment for model shot=-1 or preparation shot >=100'000
+% special treatment for model shot=-1 or preparation shot >=900'000
 begstr = '';
 if (iscell(mapping_for_tcv.expression) || isempty(strfind(mapping_for_tcv.expression,'\rtc::'))) && ...
-      ~isempty(shot) && (shot==-1 || (shot>=100000 && shot < 200000) || liuqe_version==-1 )
+      ~isempty(shot) && (shot==-1 || (shot>=900000 && shot <= 999999 ) || liuqe_version==-1 )
   % requires FBTE
   liuqe_version_eff = -1;
   begstr = 'tcv_eq( "';
@@ -3613,14 +3613,14 @@ catch
     gdat_data.error_bar_raw = gdat_data.error_bar;
   end
 
-  if (nverbose>=1) && shot<100000
+  if (nverbose>=1) && shot<900000
     warning('Problems with \results::thomson:times')
     warning(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
   end
   return
 end
 if isempty(time) || ischar(time)
-  if (nverbose>=1) && shot<100000
+  if (nverbose>=1) && shot<900000
     warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
     disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
   end
diff --git a/matlab/TCV/loadTCVdata.m b/matlab/TCV/loadTCVdata.m
deleted file mode 100644
index 49f83dcfe2d9e56ae0987a296c312b64e0808a23..0000000000000000000000000000000000000000
--- a/matlab/TCV/loadTCVdata.m
+++ /dev/null
@@ -1,1424 +0,0 @@
- function [trace,error_status,varargout]=loadTCVdata(shot,data_type,varargin)
-%
-% Added option to load shot=-1 or >=100000
-% Added option shot=-9 to list keywords
-%
-% list of data_type currently available (when [_2,_3] is added, means can add _i to get Liuqe i):
-% if -1 is added, can also get it from FBTE with shot=-1, >=100000 or liuqe_version='_-1' (to get model file)
-%
-% 'Ip'[_2,_3]   =  current
-% 'B0'[_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
-% 'j_tor'[_2,_3] =  J_TOR vs (R,Z,time)
-% 'neint' =  line-integrated electron density [m/m^3]
-% 'nel' =  line-averaged electron density [1/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'[_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
-% '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)
-% 'neft:4' =  ne fitted from data on rho mesh (from proffit.local_time:neft_abs from trial:trial_indx=4)
-% '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)
-% 'teftav' =  te fitted from averaged over time data on rho mesh (from proffit.avg_time:teft)
-% 'vtor'= Toroidal rotation of C raw profile on (rho,t). ADD error bars in .std, as well as time and rho error bars
-% 'vtorfit'= Toroidal rotation of C fitted profile on (rho,t)
-% 'vpol'= Poloidal rotation of C raw profile on (rho,t). ADD error bars in .std, as well as time and rho error bars
-% 'vpolfit'= Poloidal rotation of C fitted profile on (rho,t)
-% 'Ti' (or Tc), 'Tifit', 'ni' (or nC), 'nifit', 'zeffcxrs', 'zeffcxrsfit': similar to 'vtor' from CXRS
-% '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
-% 'IOH'  =  current in ohmic transformer coils IOH_1
-% 'vloop'  =  loop voltage
-% 'pgyro' =  ECH power for each gyro(1:9) and total (10)
-%
-% 'xx_2 or xx_3' for Liuqe2 or 3: same as above for xx related to equilibrium
-% 'xx_-1 or xx_-1' for FBTE values (model or shot=-1 or shot>=100000): 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}:  [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
-%                       (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:
-%
-% 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
-% trace.units:  units of data
-%
-% 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_status]=loadTCVdata(shot,'Ip',1);
-%         [sxr,error_status,R]=loadTCVdata(shot,'sxR',1);
-%
-
-varargout{1}=cell(1,1);
-error_status=1;
-
-% To allow multiple ways of writing a specific keyword, use data_type_eff within this routine
-if exist('data_type') && ~isempty(data_type)
-  data_type_eff=data_type;
-else
-  data_type_eff=' ';
-end
-i_23=0;
-% LIUQE tree
-begstr = '\results::';
-endstr = '';
-liuqe_version = 1;
-if length(data_type_eff)>2
-  if strcmp(data_type_eff(end-1:end),'_2') | strcmp(data_type_eff(end-1:end),'_3')
-    i_23=2;
-    endstr=data_type_eff(end-1:end);
-    liuqe_version = str2num(data_type_eff(end:end));
-  elseif strcmp(upper(data_type_eff(end-2:end)),'_-1')
-    i_23=3;
-    begstr = 'tcv_eq( "';
-    endstr = '", "FBTE" )';
-    liuqe_version = -1;
-  end
-end
-if shot==-1 || (shot>=100000 && shot<200000)
-  % requires FBTE
-  liuqe_version = -1;
-  begstr = 'tcv_eq( "';
-  endstr = '", "FBTE" )';
-end
-
-% use keyword without eventual _2 or _3 extension to check for multiple possibilities
-% Also remove ":4" for trial_indx specification
-jj=strfind(data_type_eff,':');
-trialindx=[];
-if ~isempty(jj)
-  ii=strmatch(data_type_eff(1:jj-1),[{'teft'},{'neft'},{'teftav'},{'neftav'}]);
-  if ~isempty(ii)
-    trialindx=str2num(data_type_eff(jj+1:end));
-    data_type_eff_noext=[data_type_eff(1:jj-1) ':trial'];
-  else
-    data_type_eff_noext=data_type_eff(1:end-i_23);
-  end
-else
-  data_type_eff_noext=data_type_eff(1:end-i_23);
-end
-
-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,[{'b0'} {'B_0'} {'Bgeom'} {'BGEOM'}],'exact'))
-  data_type_eff_noext='B0';
-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,[{'Te_edge'} {'TE_edge'}],'exact'))
-  data_type_eff_noext='te_edge';
-end
-if ~isempty(strmatch(data_type_eff_noext,[{'Ne_edge'} {'NE_edge'}],'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,[{'Terhozshift'}],'exact'))
-  data_type_eff_noext='terhozshift';
-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,[{'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
-if ~isempty(strmatch(data_type_eff_noext,[{'elongation'} {'elon'}],'exact'))
-  data_type_eff_noext='kappa';
-end
-if ~isempty(strmatch(lower(data_type_eff_noext),[{'j_tor'} {'jtor'} {'\results::j_tor'}],'exact'))
-  data_type_eff_noext='jtor';
-  data_type_eff = ['jtor' endstr];
-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,[{'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
-if ~isempty(strmatch(data_type_eff_noext,[{'vtor'} {'v_tor'} {'vi'} {'vc'} {'v_i'} {'v_c'} {'VTOR'} {'V_TOR'} {'VI'} {'VC'} {'V_I'} {'V_C'}],'exact'))
-  data_type_eff_noext='vi_tor';
-end
-if ~isempty(strmatch(data_type_eff_noext,[{'vtorfit'} {'vtorft'} {'v_torfit'} {'vifit'} {'vcfit'} {'v_ifit'} {'v_cfit'} {'VTORfit'} {'V_TORfit'} {'VIfit'} {'VCfit'} {'V_Ifit'} {'V_Cfit'}],'exact'))
-  data_type_eff_noext='vi_torfit';
-end
-if ~isempty(strmatch(data_type_eff_noext,[{'vpol'} {'v_pol'} {'VPOL'} {'V_POL'}],'exact'))
-  data_type_eff_noext='vi_pol';
-end
-if ~isempty(strmatch(data_type_eff_noext,[{'vpolfit'} {'vpolft'} {'v_polfit'} {'VPOLfit'} {'V_POLfit'}],'exact'))
-  data_type_eff_noext='vi_polfit';
-end
-if ~isempty(strmatch(data_type_eff_noext,[{'Ti'} {'ti'} {'TC'} {'tc'} {'T_C'} {'t_c'}],'exact'))
-  data_type_eff_noext='Ti';
-end
-if ~isempty(strmatch(data_type_eff_noext,[{'Tifit'} {'tifit'} {'TCfit'} {'tcfit'} {'T_Cfit'} {'t_cfit'}],'exact'))
-  data_type_eff_noext='Tifit';
-end
-if ~isempty(strmatch(data_type_eff_noext,[{'Ni'} {'ni'} {'NC'} {'nc'} {'N_C'} {'n_c'}],'exact'))
-  data_type_eff_noext='ni';
-end
-if ~isempty(strmatch(data_type_eff_noext,[{'Nifit'} {'nifit'} {'NCfit'} {'ncfit'} {'N_Cfit'} {'n_cfit'}],'exact'))
-  data_type_eff_noext='nifit';
-end
-if ~isempty(strmatch(data_type_eff_noext,[{'Zeffcxrs'} {'zeffcxrs'} {'Z_effcxrs'} {'z_effcxrs'} {'Zeff_cxrs'} {'zeff_cxrs'} {'Z_eff_cxrs'} {'z_eff_cxrs'}],'exact'))
-  data_type_eff_noext='zeffcxrs';
-end
-if ~isempty(strmatch(data_type_eff_noext,[{'Zeffcxrsfit'} {'zeffcxrsfit'} {'Z_effcxrsfit'} {'z_effcxrsfit'} {'Zeff_cxrsfit'} {'zeff_cxrsfit'} {'Z_eff_cxrsfit'} {'z_eff_cxrsfit'}],'exact'))
-  data_type_eff_noext='zeffcxrsfit';
-end
-
-% some defaults
-
-% nodes which have _2 and _3 equivalence, related to simpletdi case
-liuqe23=[{'\results::area'} {'\results::beta_pol'} {'\results::beta_tor'} ...
-	 {'\results::delta_95'} {'\results::delta_95_bot'} {'\results::delta_95_top'} ...
-	 {'\results::delta_edge'} {'\results::delta_ed_bot'} {'\results::delta_ed_top'} ...
-	 {'\results::diamag'} {'\results::i_p'} {'\results::j_tor'} ...
-	 {'\results::kappa_95'} {'\results::kappa_edge'} {'\results::kappa_zero'} ...
-	 {'\results::lambda'} {'\results::l_i'} {'\results::psi_axis'} {'\results::psi_values'} ...
-	 {'\results::q_95'} {'\results::q_edge'} {'\results::q_zero'} {'\results::rms_error'} ...
-	 {'\results::z_axis'} {'\results::r_axis'} {'\results::q_psi'} ...
-         {'\results::total_energy'} {'\results::tor_flux'} {'\results::volume'} ...
-         {'\results::r_contour'} {'\results::z_contour'} ...
-         {'\results::thomson:psiscatvol'} {'\results::thomson:psi_max'}];
-
-% nodes which have FBTE equivalence, related to simpletdi case
-liuqeFBTE=[{'\results::i_p'} {'\results::z_axis'} {'\results::r_axis'} {'\results::q_psi'} ...
-	   {'\results::beta_tor'} {'\results::beta_pol'} {'\results::q_95'} {'\results::l_i'}  {'\results::delta_95'} ...
-	   {'\results::kappa_95'} {'\results::r_contour'} {'\results::z_contour'} {'\results::psi_axis'}];
-
-% keywords which do not have FBTE equivalence and are returned empty
-noFBTE=[{'ne'} {'te'} {'nerho'} {'terho'} {'ne_edge'} {'te_edge'} {'nerho_edge'} {'terho_edge'} {'nerhozshift'} {'terhozshift'} {'profnerho'} {'profterho'} ...
-	{'neft'} {'neft:trial'} {'teft:trial'} {'neftav:trial'} {'teftav:trial'} {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'} {'neint'} {'nel'} ...
-        {'vi_tor'} {'vi_torfit'} {'vi_pol'} {'vi_polfit'} {'Ti'} {'Tifit'} {'ni'} {'nifit'} {'zeffcxrs'} {'zeffcxrsfit'}];
-
-% all keywords and corresponding case to run below
-TCVkeywrdall=[{'Ip'} {'B0'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'rhovol'} {'qrho'} {'q95'} {'kappa'} ...
-      {'delta'} {'deltatop'} {'deltabot'} {'neint'} {'nel'} ...
-      {'ne'} {'te'} {'nerho'} {'terho'}  {'ne_edge'} {'te_edge'} {'nerho_edge'} {'terho_edge'} {'nerhozshift'} {'terhozshift'} {'profnerho'} {'profterho'} ...
-      {'neft'} {'teft'} {'neftav'} {'teftav'} {'neft:trial'} {'teft:trial'} {'neftav:trial'} {'teftav:trial'}  ...
-      {'sxr'} {'sxR'} {'ece'} {'MPX'} {'IOH'} {'vloop'} {'pgyro'} {'jtor'} {'vi_tor'} {'vi_torfit'} {'vi_pol'} {'vi_polfit'} {'Ti'} {'Tifit'} {'ni'} {'nifit'} {'zeffcxrs'} {'zeffcxrsfit'}];
-
-TCVsig.iip=strmatch('Ip',TCVkeywrdall,'exact');
-TCVsig.iB0=strmatch('B0',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');
-TCVsig.idelta=strmatch('delta',TCVkeywrdall,'exact');
-TCVsig.ideltatop=strmatch('deltatop',TCVkeywrdall,'exact');
-TCVsig.ideltabot=strmatch('deltabot',TCVkeywrdall,'exact');
-TCVsig.ineint=strmatch('neint',TCVkeywrdall,'exact');
-TCVsig.inel=strmatch('nel',TCVkeywrdall,'exact');
-TCVsig.ine=strmatch('ne',TCVkeywrdall,'exact');
-TCVsig.ite=strmatch('te',TCVkeywrdall,'exact');
-TCVsig.ine_edge=strmatch('ne_edge',TCVkeywrdall,'exact');
-TCVsig.ite_edge=strmatch('te_edge',TCVkeywrdall,'exact');
-TCVsig.inerho=strmatch('nerho',TCVkeywrdall,'exact');
-TCVsig.iterho=strmatch('terho',TCVkeywrdall,'exact');
-TCVsig.inerho_edge=strmatch('nerho_edge',TCVkeywrdall,'exact');
-TCVsig.iterho_edge=strmatch('terho_edge',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');
-TCVsig.ineft_trial=strmatch('neft:trial',TCVkeywrdall);
-TCVsig.iteft=strmatch('teft',TCVkeywrdall,'exact');
-TCVsig.iteft_trial=strmatch('teft:trial',TCVkeywrdall);
-TCVsig.ineftav=strmatch('neftav',TCVkeywrdall,'exact');
-TCVsig.ineftav_trial=strmatch('neftav:trial',TCVkeywrdall);
-TCVsig.iteftav=strmatch('teftav',TCVkeywrdall,'exact');
-TCVsig.iteftav_trial=strmatch('teftav:trial',TCVkeywrdall);
-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');
-TCVsig.ipgyro=strmatch('pgyro',TCVkeywrdall,'exact');
-TCVsig.ijtor=strmatch('jtor',TCVkeywrdall,'exact');
-TCVsig.ivi_tor=strmatch('vi_tor',TCVkeywrdall,'exact');
-TCVsig.ivi_torfit=strmatch('vi_torfit',TCVkeywrdall,'exact');
-TCVsig.ivi_pol=strmatch('vi_pol',TCVkeywrdall,'exact');
-TCVsig.ivi_polfit=strmatch('vi_polfit',TCVkeywrdall,'exact');
-TCVsig.iTi=strmatch('Ti',TCVkeywrdall,'exact');
-TCVsig.iTifit=strmatch('Tifit',TCVkeywrdall,'exact');
-TCVsig.ini=strmatch('ni',TCVkeywrdall,'exact');
-TCVsig.inifit=strmatch('nifit',TCVkeywrdall,'exact');
-TCVsig.izeffcxrs=strmatch('zeffcxrs',TCVkeywrdall,'exact');
-TCVsig.izeffcxrsfit=strmatch('zeffcxrsfit',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.iB0)=TCVkeywrdall(TCVsig.iB0); % through iphi
-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.ine_edge)=TCVkeywrdall(TCVsig.ine_edge); % special as dimensions from other nodes
-TCVkeywrdcase(TCVsig.ite_edge)=TCVkeywrdall(TCVsig.ite_edge); % idem
-TCVkeywrdcase(TCVsig.inerho)=TCVkeywrdall(TCVsig.inerho); % idem
-TCVkeywrdcase(TCVsig.inerho_edge)=TCVkeywrdall(TCVsig.inerho_edge); % idem
-TCVkeywrdcase(TCVsig.iteft_trial)=TCVkeywrdall(TCVsig.iteft_trial); % special to extract trial_indx if needed
-TCVkeywrdcase(TCVsig.ineft_trial)=TCVkeywrdall(TCVsig.ineft_trial); % special to extract trial_indx if needed
-TCVkeywrdcase(TCVsig.iteftav_trial)=TCVkeywrdall(TCVsig.iteftav_trial); % special to extract trial_indx if needed
-TCVkeywrdcase(TCVsig.ineftav_trial)=TCVkeywrdall(TCVsig.ineftav_trial); % special to extract trial_indx if needed
-TCVkeywrdcase(TCVsig.iterho)=TCVkeywrdall(TCVsig.iterho); % idem
-TCVkeywrdcase(TCVsig.iterho_edge)=TCVkeywrdall(TCVsig.iterho_edge); % 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);
-TCVkeywrdcase(TCVsig.iMPX)=TCVkeywrdall(TCVsig.iMPX);
-TCVkeywrdcase(TCVsig.iIOH)=TCVkeywrdall(TCVsig.iIOH);
-TCVkeywrdcase(TCVsig.ivloop)=TCVkeywrdall(TCVsig.ivloop);
-TCVkeywrdcase(TCVsig.ipgyro)=TCVkeywrdall(TCVsig.ipgyro);
-TCVkeywrdcase(TCVsig.ijtor)=TCVkeywrdall(TCVsig.ijtor);
-TCVkeywrdcase(TCVsig.ivi_tor)=TCVkeywrdall(TCVsig.ivi_tor);
-TCVkeywrdcase(TCVsig.ivi_torfit)=TCVkeywrdall(TCVsig.ivi_torfit);
-TCVkeywrdcase(TCVsig.ivi_pol)=TCVkeywrdall(TCVsig.ivi_pol);
-TCVkeywrdcase(TCVsig.ivi_polfit)=TCVkeywrdall(TCVsig.ivi_polfit);
-TCVkeywrdcase(TCVsig.iTi)=TCVkeywrdall(TCVsig.iTi);
-TCVkeywrdcase(TCVsig.iTifit)=TCVkeywrdall(TCVsig.iTifit);
-TCVkeywrdcase(TCVsig.ini)=TCVkeywrdall(TCVsig.ini);
-TCVkeywrdcase(TCVsig.inifit)=TCVkeywrdall(TCVsig.inifit);
-TCVkeywrdcase(TCVsig.izeffcxrs)=TCVkeywrdall(TCVsig.izeffcxrs);
-TCVkeywrdcase(TCVsig.izeffcxrsfit)=TCVkeywrdall(TCVsig.izeffcxrsfit);
-
-% 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.iq95)={'\results::q_95'};
-TCVsiglocation(TCVsig.ikappa)={'\results::kappa_edge'};
-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.inel)={'\results::fir:n_average'};
-%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;
-TCVsiglocation(TCVsig.iteftav)={'\results::proffit.avg_time:teft'}; TCVsigtimeindx(TCVsig.iteftav)=2;
-TCVsiglocation(TCVsig.ineft_trial)=TCVsiglocation(TCVsig.ineft); TCVsigtimeindx(TCVsig.ineft_trial)=2;
-TCVsiglocation(TCVsig.iteft_trial)=TCVsiglocation(TCVsig.iteft); TCVsigtimeindx(TCVsig.iteft_trial)=2;
-TCVsiglocation(TCVsig.ineftav_trial)=TCVsiglocation(TCVsig.ineftav); TCVsigtimeindx(TCVsig.ineftav_trial)=2;
-TCVsiglocation(TCVsig.iteftav_trial)=TCVsiglocation(TCVsig.iteftav); TCVsigtimeindx(TCVsig.iteftav_trial)=2;
-
-if shot==-9
-  clear trace
-  for i=1:length(TCVkeywrdall)
-    fieldname_eff = lower(TCVkeywrdall{i});
-    keyword_eff = TCVkeywrdall{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 = TCVsiglocation{TCVsig.iip};
-  % trace.data=[];
-  return
-end
-
-% initialize order of substructures and allows just a "return" if data empty
-trace.data=[];
-trace.x=[];
-trace.t=[];
-trace.dim=[];
-trace.dimunits=[];
-trace.units=[];
-iprintwarn=0;
-% 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 iprintwarn & isempty(index)
-    disp('********************')
-    disp('trace not yet registered.')
-    disp('If standard data, ask 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('********************')
-  elseif isempty(index)
-    % 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_noext,TCVkeywrdall,'exact');
-  if isempty(index)
-    disp(' ')
-    disp('********************')
-    disp(['no such keyword: ' data_type_eff(1:end-i_23)])
-    disp(' ')
-    disp('Available keywords:')
-    TCVkeywrdall(:)
-    disp('********************')
-    return
-  end
-  if liuqe_version==-1 && ~isempty(strmatch(TCVkeywrdall{index},noFBTE,'exact'))
-    disp(['node ' TCVkeywrdall{index} ' not defined within FBTE'])
-    return
-  end
-end
-if iprintwarn
-  disp(['loading' ' ' data_type_eff_noext ' from TCV shot #' num2str(shot)]);
-  disp(['case ' TCVkeywrdcase{index}])
-end
-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'])
-  disp('********')
-end
-
-status=ones(1,100);
-zmag=cell(0,0);
-nargineff=nargin;
-if iprintwarn
-  disp('TCVsiglocation{index} in loadTCVdata')
-  TCVsiglocation{index}
-end
-
-switch TCVkeywrdcase{index}
-
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case 'simpletdi'
-
-    %  load TCV other data
-    if liuqe_version==-1
-      mdsopen( 'pcs', shot); %synthetic shot generated with FBTE and MGAMS.
-    else
-      mdsopen(shot);
-      % test if node exists
-      error_status=1;
-      ij=findstr(TCVsiglocation{index},'[');
-      if isempty(ij); ij=length(TCVsiglocation{index})+1; end
-      if eval(['~mdsdata(''node_exists("\' TCVsiglocation{index}(1:ij-1) '")'')'])
-	disp(['node ' TCVsiglocation{index}(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)])
-% $$$       return
-      end
-    end
-    if strcmp(TCVsiglocation{index}(1:9),'\psitbx::')
-      begstr = 'tcv_psitbx("';
-      nodenameeff=[begstr TCVsiglocation{index}(10:end) endstr];
-    elseif length(TCVsiglocation{index})>16 && strcmp(TCVsiglocation{index}(1:16),'\results::psitbx')
-      begstr = 'tcv_psitbx("';
-      nodenameeff=[begstr TCVsiglocation{index}(18:end) endstr];
-    elseif  strcmp(TCVsiglocation{index}(1:10),'\results::')
-      nodenameeff=[begstr TCVsiglocation{index}(11:end) endstr];
-    else
-      nodenameeff=TCVsiglocation{index};
-      disp(['should not have gone through here, mention this example to O. Sauter']);
-    end
-    tracetdi=tdi(nodenameeff);
-    mdsclose;
-    if isempty(tracetdi.data) || (~iscell(tracetdi.data) && (length(tracetdi.data)<=1) && isnan(tracetdi.data))
-      disp(['node ' nodenameeff ' is empty for shot = ' num2str(shot)])
-      return
-    end
-    trace.data=tracetdi.data;
-
-    if length(tracetdi.dim)==0
-      % scalar
-      trace.x=[];
-      trace.t=[];      
-      elseif 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
-      % if dim=41 assumes is x (usual for TCV)
-      if length(tracetdi.dim{1})==41
-        ix=1;
-        trace.t=tracetdi.dim{2};
-      elseif length(tracetdi.dim{2})==41
-        ix=2;
-        trace.t=tracetdi.dim{1};
-      elseif length(tracetdi.dim)==2
-        % if one has 2D or more and the other 1D, assume rho is 2D or more and time is 1D
-        if min(size(tracetdi.dim{1}))==1 && (min(size(tracetdi.dim{2}))>1 || length(size(tracetdi.dim{2}))>2)
-          ix=2;
-          trace.t=tracetdi.dim{1};
-          disp(['assumes time array is 1D array, so x from dim{' num2str(ix) '}'])
-        elseif min(size(tracetdi.dim{2}))==1 && (min(size(tracetdi.dim{1}))>1 || length(size(tracetdi.dim{1}))>2)
-          ix=1; 
-          trace.t=tracetdi.dim{2};
-          disp(['assumes time array is 1D array, so x from dim{' num2str(ix) '}'])
-        elseif tracetdi.dim{1}(end)<1.1 && tracetdi.dim{2}(end)>1.1
-          ix=1; 
-          trace.t=tracetdi.dim{2};
-          disp(['assumes x array is 1D array with max <1.1'])          
-        elseif tracetdi.dim{1}(end)>1.1 && tracetdi.dim{2}(end)<1.1
-          ix=2; 
-          trace.t=tracetdi.dim{1};
-          disp(['assumes x array is 1D array with max <1.1'])          
-        elseif length(size(tracetdi.dim{1}))==2 && min(size(tracetdi.dim{1}))==1 && length(size(tracetdi.dim{2}))==2 && min(size(tracetdi.dim{2}))==1
-          % if both are 1D arrays, 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
-      end
-    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
-      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=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;
-    end
-    % 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];
-    error_status=0;
-
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case {'ne','te'}
-    % ne or Te from Thomson data on raw z mesh vs (z,t)
-    mdsopen(shot);
-    if strcmp(TCVkeywrdcase{index},'ne')
-      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');
-      trace_abs=tdi('\results::thomson:fir_thom_rat');
-    end
-    trace.data=tracetdi.data'; % Thomson data as (t,z)
-    trace.std=tracestd.data';
-    trace.name=[num2str(shot) ';' nodenameeff];
-    % add correct dimensions
-    try
-      time=mdsdata('\results::thomson:times');
-    catch
-      warning('Problems with \results::thomson:times')
-      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}])
-      return
-    end
-    if isempty(time) || ischar(time)
-      thomsontimes=time
-      warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
-      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}])
-      return
-    end
-    z=mdsdata('\diagz::thomson_set_up:vertical_pos');
-    trace.dim=[{z};{time}];
-    trace.dimunits=[{'Z [m]'} ; {'time [s]'}];
-    trace.x=z;
-    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('mdsip');
-
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case {'ne_edge','te_edge'}
-    % ne or Te from Thomson.edge data on raw z mesh vs (z,t)
-    mdsopen(shot);
-    if strcmp(TCVkeywrdcase{index},'ne_edge')
-      nodenameeff='\results::thomson.edge:ne';
-      tracetdi=tdi(nodenameeff);
-      tracestd=tdi('\results::thomson.edge:ne:error_bar');
-    else
-      nodenameeff='\results::thomson.edge:te';
-      tracetdi=tdi(nodenameeff);
-      tracestd=tdi('\results::thomson.edge:te:error_bar');
-      trace_abs=tdi('\results::thomson.edge:fir_thom_rat');
-    end
-    trace.data=tracetdi.data'; % Thomson.Edge data as (t,z)
-    trace.std=tracestd.data';
-    trace.name=[num2str(shot) ';' nodenameeff];
-    % add correct dimensions
-    try
-      time=mdsdata('\results::thomson:times');
-    catch
-      warning('Problems with \results::thomson:times')
-      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}])
-      return
-    end
-    if isempty(time) || ischar(time)
-      thomsontimes=time
-      warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
-      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}])
-      return
-    end
-    z=mdsdata('\diagz::thomson_set_up.edge:vertical_pos');
-    trace.dim=[{z};{time}];
-    trace.dimunits=[{'Z [m]'} ; {'time [s]'}];
-    trace.x=z;
-    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('mdsip');
-
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case {'nerho','terho'}
-    % ne or Te from Thomson data on rho=sqrt(psi_normalised) mesh: (rho,t)
-    mdsopen(shot);
-    try
-      time=mdsdata('\results::thomson:times');
-    catch
-      warning('Problems with \results::thomson:times')
-      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}])
-      return
-    end
-    if isempty(time) || ischar(time)
-      thomsontimes=time
-      warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
-      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}])
-      return
-    end
-    if strcmp(TCVkeywrdcase{index},'nerho')
-      nodenameeff='\results::thomson:ne';
-      tracetdi=tdi(nodenameeff);
-      if isempty(tracetdi.data)
-        ishot=mdsopen(shot)
-        tracetdi=tdi(nodenameeff)
-      end
-      nodenameeff='\results::thomson:ne; error_bar ; fir_thom_rat; (ne,std)*fir_thom_rat';
-      tracestd=tdi('\results::thomson:ne:error_bar');
-      if shot>=23801
-        tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
-        if isempty(tracefirrat.data) || ischar(tracefirrat.data)
-          disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty')
-        end
-      else
-        tracefirrat=tdi('\results::thomson:fir_thom_rat');
-        if isempty(tracefirrat.data) || ischar(tracefirrat.data)
-          disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty')
-          tracefirrat.dim{1}=[];
-        else
-          tracefirrat.dim{1}=time;
-        end
-      end
-      if ~isempty(tracefirrat.data) || ischar(tracefirrat.data)
-        tracefirrat_data=NaN*ones(size(tracetdi.dim{1}));
-        itim=iround(time,tracefirrat.dim{1});
-        tracefirrat_data(itim)=tracefirrat.data;
-      else
-        tracefirrat_data=NaN;
-      end
-    else
-      nodenameeff='\results::thomson:te';
-      tracetdi=tdi(nodenameeff);
-      nodenameeff='\results::thomson:te; error_bar';
-      tracestd=tdi('\results::thomson:te:error_bar');
-    end
-    trace.data=tracetdi.data'; % Thomson data as (t,z)
-    trace.std=tracestd.data';
-    if strcmp(TCVkeywrdcase{index},'nerho')
-      trace.firrat=tracefirrat_data;
-      trace.data_abs=trace.data*diag(tracefirrat_data);
-      trace.std_abs=trace.std*diag(tracefirrat_data);
-    end
-    trace.name=[num2str(shot) ';' nodenameeff];
-    % add correct dimensions
-    % construct rho mesh
-    psi_max=tdi(['\results::thomson:psi_max' endstr]);
-    psiscatvol=tdi(['\results::thomson:psiscatvol' endstr]);
-    if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data)
-      for ir=1:length(psiscatvol.dim{2})
-        rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))';
-      end
-    else
-      rho=NaN;
-    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 {'nerho_edge','terho_edge'}
-    % ne or Te from Thomson.Edge data on rho=sqrt(psi_normalised) mesh: (rho,t)
-    mdsopen(shot);
-    try
-      time=mdsdata('\results::thomson:times');
-    catch
-      warning('Problems with \results::thomson:times')
-      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}])
-      return
-    end
-    if isempty(time) || ischar(time)
-      thomsontimes=time
-      warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
-      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}])
-      return
-    end
-    if strcmp(TCVkeywrdcase{index},'nerho_edge')
-      nodenameeff='\results::thomson.edge:ne';
-      tracetdi=tdi(nodenameeff);
-      if isempty(tracetdi.data) || ischar(tracetdi.data)
-        ishot=mdsopen(shot)
-        tracetdi=tdi(nodenameeff)
-      end
-      nodenameeff='\results::thomson.edge:ne; error_bar ; fir_thom_rat; (ne,std)*fir_thom_rat';
-      tracestd=tdi('\results::thomson.edge:ne:error_bar');
-      if shot>=23801
-        tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
-        if isempty(tracefirrat.data) || ischar(tracefirrat.data)
-          disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty')
-        end
-      else
-        tracefirrat=tdi('\results::thomson.edge:fir_thom_rat');
-        tracefirrat.dim{1}=time;
-      end
-      if ~isempty(tracetdi.dim)
-        tracefirrat_data=NaN*ones(size(tracetdi.dim{1}));
-      else
-        tracefirrat_data=[];
-      end
-      if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data)
-        itim=iround(time,tracefirrat.dim{1});
-        tracefirrat_data(itim)=tracefirrat.data;
-      end
-    else
-      nodenameeff='\results::thomson.edge:te';
-      tracetdi=tdi(nodenameeff);
-      nodenameeff='\results::thomson.edge:te; error_bar';
-      tracestd=tdi('\results::thomson.edge:te:error_bar');
-    end
-    if ~ischar(tracetdi.data); trace.data=tracetdi.data'; end % Thomson.Edge data as (t,z)
-    trace.std=tracestd.data';
-    if strcmp(TCVkeywrdcase{index},'nerho')
-      trace.firrat=tracefirrat_data;
-      trace.data_abs=trace.data*diag(tracefirrat_data);
-      trace.std_abs=trace.std*diag(tracefirrat_data);
-    end
-    trace.name=[num2str(shot) ';' nodenameeff];
-    % add correct dimensions
-    % construct rho mesh
-    psi_max=tdi(['\results::thomson.edge:psi_max' endstr]);
-    psiscatvol=tdi(['\results::thomson.edge:psiscatvol' endstr]);
-    if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data)
-      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
-    else
-      trace.dim={};
-      trace.dimunits={};
-      trace.x=[];
-      trace.t=[];
-      trace.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>=3 & ~isempty(varargin{1})
-      zshift=varargin{1};
-    else
-      zshift=0.;
-    end
-    mdsopen(shot);
-    try
-      time=mdsdata('\results::thomson:times');
-    catch
-      warning('Problems with \results::thomson:times')
-      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}])
-      return
-    end
-    if isempty(time) || ischar(time)
-      thomsontimes=time
-      warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
-      disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' TCVkeywrdcase{index}])
-      return
-    end
-    if strcmp(TCVkeywrdcase{index},'nerhozshift')
-      nodenameeff='\results::thomson:ne';
-      tracetdi=tdi(nodenameeff);
-      tracestd=tdi('\results::thomson:ne:error_bar');
-      if shot>=23801
-        tracefirrat=tdi('\results::thomson.profiles.auto:fir_thom_rat'); %time base not same!!
-        if isempty(tracefirrat.data) || ischar(tracefirrat.data)
-          disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty')
-        end
-      else
-        tracefirrat=tdi('\results::thomson:fir_thom_rat');
-        if isempty(tracefirrat.data) || ischar(tracefirrat.data)
-          disp('problem with \results::thomson.profiles.auto:fir_thom_rat: empty')
-          tracefirrat.dim{1}=[];
-        else
-          tracefirrat.dim{1}=time;
-        end
-      end
-      if ~isempty(tracefirrat.data) || ischar(tracefirrat.data)
-        tracefirrat_data=NaN*ones(size(tracetdi.dim{1}));
-        itim=iround(time,tracefirrat.dim{1});
-        tracefirrat_data(itim)=tracefirrat.data;
-      else
-        tracefirrat_data=NaN;
-      end
-    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';
-    if strcmp(TCVkeywrdcase{index},'nerhozshift')
-      trace.firrat=tracefirrat_data;
-      trace.data_abs=trace.data*diag(tracefirrat_data);
-      trace.std_abs=trace.std*diag(tracefirrat_data);
-    end
-    trace.name=[num2str(shot) ';' nodenameeff];
-    % add correct dimensions
-    % construct rho mesh
-    if strcmp(endstr,'_-1')
-      error(['in ' TCVkeywrdcase{index} ' endstr should not be ' endstr]);
-    end
-    psi_max=tdi(['\results::thomson:psi_max' endstr]);
-    psiscatvol=tdi(['\results::thomson:psiscatvol' endstr]);
-    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
-    if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data)
-      for ir=1:length(psiscatvol.dim{2})
-        rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))';
-      end
-    else
-      rho=NaN;
-    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);
-    error_status=1;
-    if strcmp(TCVkeywrdcase{index},'profnerho')
-      nodenameeff=['\results::THOMSON.PROFILES.AUTO:ne'];
-      avers=tdi('\results::THOMSON.PROFILES.AUTO:ne:version_num');
-    end
-    if strcmp(TCVkeywrdcase{index},'profterho')
-      nodenameeff=['\results::THOMSON.PROFILES.AUTO:te'];
-      avers=tdi('\results::THOMSON.PROFILES.AUTO:te:version_num');
-    end
-    if avers.data==0
-     % may be because nodes not yet filled in, so call once a node
-     ab=tdi(nodenameeff);
-     avers=tdi([nodenameeff ':version_num']);
-    end
-    if avers.data>0
-      tracetdi=tdi(nodenameeff);
-      if avers.data < 2.99
-        % for earlier version the bug made it to have logically (rho,t)
-        if ~isempty(tracetdi.dim) && ~ischar(tracetdi.data)
-          trace.x=tracetdi.dim{1};
-          trace.t=tracetdi.dim{2};
-          error_status=0;
-        else
-         error_status=2;
-          trace.x=[];
-          trace.t=[];
-        end
-      else
-        trace.data=tracetdi.data'; % error in dimensions for autofits
-        if ~isempty(tracetdi.dim) && ~ischar(tracetdi.data)
-          disp('assumes dim{2} for x in THOMSON.PROFILES.AUTO')
-          trace.x=tracetdi.dim{2};
-          trace.t=tracetdi.dim{1};
-          error_status=0;
-        else
-          trace.x=[];
-          trace.t=[];
-         error_status=2;
-        end
-      end
-    else
-      tracetdi=avers;
-      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;
-
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case {'B0'}
-    % B0 at R0=0.88
-    if liuqe_version==-1
-      % mdsopen( 'pcs', shot); %synthetic shot generated with FBTE and MGAMS.
-      mdsopen(shot)
-      nodenameeff = 'tcv_eq("BZERO","FBTE")';
-      tracetdi=tdi(nodenameeff);
-      trace.data = tracetdi.data;
-      trace.t = tracetdi.dim{1};
-      
-    else
-      mdsopen(shot);
-      nodenameeff=['\magnetics::iphi'];
-      tracetdi=tdi(nodenameeff);
-      R0EXP=0.88;
-      trace.data=192.E-07 * 0.996 *tracetdi.data/R0EXP;
-      trace.t=tracetdi.dim{1};
-    end
-    trace.x=[];
-%keyboard
-    trace.dim=[{trace.t}];
-    trace.dimunits=[{'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) ';' nodenameeff];
-    mdsclose;
-
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case {'qrho'}
-    % q profile on psi from liuqe
-    if liuqe_version==-1
-      mdsopen( 'pcs', shot); %synthetic shot generated with FBTE and MGAMS.
-    else
-      mdsopen(shot);
-    end
-    nodenameeff=[begstr 'q_psi' endstr];
-    tracetdi=tdi(nodenameeff);
-    trace.data=tracetdi.data;
-    trace.x=sqrt(tracetdi.dim{1}/(length(tracetdi.dim{1})-1));
-    trace.t=tracetdi.dim{2};
-    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) ';' nodenameeff];
-    mdsclose;
-
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case {'vol'}
-    % vol from psitbx
-    if liuqe_version==-1
-      begstr = 'tcv_psitbx("';
-      mdsopen( 'pcs', shot); %synthetic shot generated with FBTE and MGAMS.
-      nodenameeff=[begstr 'vol' endstr];
-      tracetdi=tdi(nodenameeff);
-    else
-      mdsopen(shot);
-      nodenameeff=['\results::psitbx:vol'];
-      tracetdi=tdi(nodenameeff);
-      if i_23==2
-	['LIUQE' endstr(2:2)]
-	psi=psitbxtcv(shot,'+0',['LIUQE' endstr(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
-    end
-    trace.data=tracetdi.data;
-    if isempty(tracetdi.data) || ischar(tracetdi.data)
-      trace.x=tracetdi.dim;
-      trace.t=tracetdi.dim;
-      trace.dim=tracetdi.dim;
-    else
-      trace.x=tracetdi.dim{1};
-      trace.t=tracetdi.dim{2};
-      trace.dim=[{trace.x};{trace.t}];
-    end
-    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) ';' nodenameeff '_' num2str(liuqe_version)];
-    mdsclose;
-
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case {'rhovol'}
-    % vol from psitbx
-    if liuqe_version==-1
-      begstr = 'tcv_psitbx("';
-      mdsopen( 'pcs', shot); %synthetic shot generated with FBTE and MGAMS.
-      nodenameeff=[begstr 'vol' endstr];
-      tracetdi=tdi(nodenameeff);
-    else
-      mdsopen(shot);
-      nodenameeff=['\results::psitbx:vol'];
-      tracetdi=tdi(nodenameeff);
-      if i_23==2
-	['LIUQE' endstr(2:2)]
-	psi=psitbxtcv(shot,'+0',['LIUQE' endstr(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
-    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]'}];
-    % 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 '_' num2str(liuqe_version)];
-    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
-    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 nargin>=5 & ~isempty(varargin{3})
-          zmag=varargin{3};
-        else
-          zmag=loadTCVdata(shot,['zmag' '_' num2str(liuqe_version)]);
-        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.data=VsxrTCVradius(data,xchord,ychord)';
-	radius.t=t';
-	varargout{1}={radius};
-        trace.R=radius.data;
-      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
-    trace.data=xtomo_signal;
-    trace.x=[1:size(trace.data,1)]';
-    trace.t=t';
-    trace.dim=[{trace.x} ; {trace.t}];
-    trace.dimunits=[{'channel #'} ; {'time [s]'}];
-    trace.name=[num2str(shot) ';' 'get_xtomo_data'];
-    error_status=0;
-
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case 'ece'
-    %  load TCV ECE data
-    % Status=1 => Not Read Yet
-    mdsopen(shot);
-    if ~isempty(find(status == 1))
-      if eval(['~mdsdata(''node_exists("\\RESULTS::ECE:rho")'')'])
-        disp(['node \RESULTS::ECE:rho does not exist for shot = ' num2str(shot)])
-        return
-      end
-      if eval(['mdsdata(''getnci("\\RESULTS::ECE:rho","length")'')==0'])
-        disp(['no data for \RESULTS::ECE:rho for shot = ' num2str(shot)])
-        return
-      end
-      [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));
-    trace.data=TE_ECE(a:b,:)';
-    trace.t=T(a:b);
-    trace.x=[1:size(trace.data,1)]';
-    trace.dim=[{trace.x} ; {trace.t}];
-    trace.dimunits=[{'channel #'} ; {'time [s]'}];
-    trace.R=R(a:b,:)';
-    trace.name=[num2str(shot) ';' '\results::ece...'];
-    radius.data=trace.R;
-    radius.t=trace.t;
-    varargout{1}={radius};
-    error_status=0;
-    % 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 '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);
-      keyboard
-      signal=get_mds_mio('MPX',[t_1 t_2]);
-      mdsclose;
-      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]'};
-    trace.name=[num2str(shot) ';' 'get_mds_mio(MPX)'];
-    [xchord,ychord]=mpx_geometry;
-    varargout{1}={VsxrTCVradius(zmag.data,xchord,ychord)};
-    error_status=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]'};
-    % 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}];
-    mdsclose;
-    error_status=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]'};
-    % 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}];
-    mdsclose;
-    error_status=0;
-    
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case 'pgyro'
-    %  ECH power for each gyro(1:9) and total (10)
-    mdsopen(shot);
-    nodenameeff=[{'\results::toray.input:p_gyro'}];
-    tracetdi=tdi(nodenameeff{:});
-    trace.data=tracetdi.data;
-    trace.x=[];
-    trace.t=tracetdi.dim{1};
-    trace.dim=tracetdi.dim;
-    trace.dimunits=tracetdi.dimunits;
-    % 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];
-    mdsclose;
-    error_status=0;
-    
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case 'jtor'
-    %  \results::j_tor , 3-D so need to specify time dim index
-    mdsopen(shot);
-    nodenameeff=[{'\results::j_tor'} endstr];
-    tracetdi=tdi(nodenameeff{:});
-    trace.data=tracetdi.data;
-    trace.x=tracetdi.dim{1};
-    trace.t=tracetdi.dim{3};
-    trace.dim=tracetdi.dim;
-    trace.dimunits=tracetdi.dimunits;
-    % 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];
-    mdsclose;
-    error_status=0;
-    
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case {'neft:trial','teft:trial','neftav:trial','teftav:trial'}
-    %  trial indx
-    mdsopen(shot);
-    eval(['nodenameeff={''' TCVsiglocation{index} ':trial''};']);
-    tracetdi=tdi(nodenameeff{:});
-    if isempty(trialindx)
-      error('trialindx should not be empty, check call or ask O. Sauter');
-    end
-    trace.data=tracetdi.data(:,:,trialindx+1);
-    trace.x=tracetdi.dim{1};
-    trace.t=tracetdi.dim{2};
-    trace.dim=tracetdi.dim(1:2);
-    trace.dimunits=tracetdi.dimunits(1:2);
-    % 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{:} ' ; trialindx=' num2str(trialindx) ];
-    mdsclose;
-    error_status=0;
-    
-  %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
- case {'vi_tor', 'vi_torfit', 'vi_pol', 'vi_polfit', 'Ti', 'Tifit', 'ni', 'nifit', 'zeffcxrs', 'zeffcxrsfit'}
-    proffit = '';
-    kwd_eff = TCVkeywrdcase{index};
-    ii=strfind(kwd_eff,'fit');
-    if ~isempty(ii); 
-      proffit = ':proffit';
-      kwd_eff = kwd_eff(1:ii-1);
-    end
-    if strcmp(kwd_eff,'zeffcxrs');  kwd_eff = 'zeff'; end
-    mdsopen(shot);
-    eval(['nodenameeff=''\results::cxrs' proffit ':' kwd_eff endstr ''';']);
-    tracetdi=tdi(nodenameeff);
-    trace.data=tracetdi.data;
-    if length(tracetdi.dim)>1;
-      trace.x=tracetdi.dim{1};
-      trace.t=tracetdi.dim{2};
-    else
-      trace.x=[];
-      trace.t=[];
-    end
-    trace.dim=tracetdi.dim;
-    trace.dimunits=tracetdi.dimunits;
-    % 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];
-    % add error bars
-    eval(['nodenameeff=''\results::cxrs' proffit ':' kwd_eff ':err' endstr ''';']);
-    nodenameeff=[{'\results::cxrs:vi_tor:err'} endstr];
-    tracetdi=tdi(nodenameeff{:});
-    trace.std = tracetdi.data;
-    if length(tracetdi.dim)>2;
-      trace.std_rho = tracetdi.dim{1};
-      trace.std_t = tracetdi.dim{2};
-    else
-      trace.std_rho = [];
-      trace.std_t = [];
-    end
-    mdsclose;
-    error_status=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
diff --git a/matlab/TCV/tcv_requests_mapping.m b/matlab/TCV/tcv_requests_mapping.m
index 21c61bec63852c398bda186d50afa492957f6bbe..9909a953e12faf245943034c299fe7db3eee168c 100644
--- a/matlab/TCV/tcv_requests_mapping.m
+++ b/matlab/TCV/tcv_requests_mapping.m
@@ -446,6 +446,17 @@ switch lower(data_request)
   mapping.timedim = 1;
   mapping.gdat_timedim = 2;
   mapping.method = 'switchcase';
+ case 'tau_r'
+  mapping.timedim = 1;
+  mapping.gdat_timedim = 1;
+  mapping.method = 'expression';
+  mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''a_minor''; ' ...
+                    'gdat_tmp=gdat_tcv([],params_eff); ' ...
+                    'params_eff = gdat_data.gdat_params;params_eff.data_request=''\results::ibs:signeo''; ' ...
+                    'aa=gdat_tcv([],params_eff);ix=iround_os(aa.x,0.4); ' ...
+                    'bb=interpos(aa.t,aa.data(ix,:),gdat_tmp.t,-100);' ...
+                    'gdat_tmp.data = 4e-7*pi.*gdat_tmp.data.^2.*bb./1.22; gdat_tmp.label=''\tau_R'';gdat_tmp.gdat_request=''tau_r'';' ...
+                    'gdat_tmp.units=''s'';gdat_tmp.help=''tauR=mu0 a^2 / 1.22 etaneo(rhopol=0.4)'''];
  case 'te'
   mapping.timedim = 2;
   mapping.label = 'Te';
diff --git a/matlab/TCV/xtomo_geometry.m b/matlab/TCV/xtomo_geometry.m
deleted file mode 100644
index a809b9e65cca29d8df8eff407eaac96d7b7ebc4c..0000000000000000000000000000000000000000
--- a/matlab/TCV/xtomo_geometry.m
+++ /dev/null
@@ -1,402 +0,0 @@
-function [fans,vangle,xchord,ychord,aomega,angfact]=xtomo_geometry(i_detec,fans)
-
-% ----[anton.public]
-%
-% function 
-%
-% [fans,vangle,xchord,ychord,aomega,angfact]=xtomo_geometry(i_detec,fans);
-% 
-%	inputs: 
-%
-%  i_detec:	=2: Xtomo prototype cameras (shot# < 6768)
-%		=1: Xtomo 9-cameras	    (shot# > 682x)
-%
-%	outputs:
-%
-%  fans:    camera switch, 1=on,0=off (1x10)
-%  vangle:  angle between detect. surface normal and pos. x-axis (1x10)
-%  xchord:  two x-coordinates (2xnl) and
-%  ychord:  two y-coord. for each line (2xnl), they specify start + end points
-%  aomega:  etendue in mm^2 x steradians
-%  angfact: angular factors, inverse of relative etendue (throughput) (20x10)
-%
-%	uses:
-%		AOMEGA=etendue_n2(b1x,b1y,b1z,b2x,b2y,b2z,z01,z02,X0,cw);
-%		angular_fact_*.mat , '*'=i_detec
-%  
-%
-%---------------- M.Anton 14/3/95 -------------------------------------------
-
-disp('*----------------------------*')
-disp('|   this is xtomo_geometry   |')
-disp('*----------------------------*')
-
-global xap yap xdet ydet
-global ae da
-% ======== tokamak parameters ================================================
-
-load tcv_vesc1
-
-
-
-
-xcont=rzvin(:,1);
-ycont=rzvin(:,2);
-xmin=min(xcont);
-ymin=min(ycont);
-xmax=max(xcont);
-ymax=max(ycont);
-xedge=100;
-yedge=60;
-          
-
-
-
-% =========  detector parameters =============================================
-
-if i_detec==2
-
-	cw=1;                         % detector numbers cw=1:clockwise cw=0:ccw
-
-	if nargin<2
-		fans=[0 0 0 0 0 0 1 0 1 0];   % camera switch
-	end
-
-	vangle=[90 90 90 0 0 0 0 -90 -90 -90];
-                                % angle of detector surface normal 
-	xpos=[0 0 0 0 0 0  118.05 0 87.84 0]; 
-                                % x position of the diaphragmas in [cm]
-	ypos=[0 0 0 0 0 0 -46 0 -80.45 0];
-                                % y position of the diaphragmas in [cm]
-	ae=[0 0 0 0 0 0 -2.5 0 -0.1 0]/10; 
-				%excentricity array/diaphragm in [cm]
-	da=[0 0 0 0 0 0 10.1  0 -11.7 0]/10;   
-                                 % diaphragma-array distance in [cm]
-	da2=da;
-
-	d1=0.950;                     % detector width      in mm
-	d2=4.050;                     % detector length     in mm
-	b1=1.000;                     % aperture width      in mm
-	b2=4.000*ones(1,10);          % aperture length     in mm
-	b3x=0;                        % aperture thickness  in mm
-	b3y=0;
-
-elseif i_detec==1
-
-   cw=1;                         % detector numbers cw=1:clockwise cw=0:ccw
-
-   if nargin<2	
-   	fans=[1 1 1 1 1 1 1 1 0 1];   % camera switch
-   end
-
-   vangle=[90 90 90 0 0 0 0 -90 -90 -90];
-                                 % angle of detector surface normal 
-   xpos=[71.5 87.7 103.9 123.1  123.1 123.1 123.1 104.04 87.84 71.64]; 
-                                 % x position of the diaphragmas in [cm]
-   ypos=[80.35 80.35 80.35 48.1 1.95 -2.45 -48.6 -80.35 -80.35 -80.35];
-                                 % y position of the diaphragmas in [cm]
-   ae=[-8 0 8 5 9 -9 -5 8 0 -8]/10; %excentricity array/diaphragm [cm]
-
-
-  ae= ae + [-0.0915 0   0.1361    0.2123    0.0923   -0.0994   ...
- 		 0.0872   -0.1520  0  0.9410 ]/10;
-  ae(1)=ae(1)+0.1/10;
-  ae(3)=8/10+0.14/10;
-  ae(4)=4.9/10;
-  ae(5)=9/10+0.2/10;
-  ae(6)=ae(6)-0.2/10;
-  ae(7)=-4.9/10;
-  ae(10)=-7.1/10;
-
-   da=   [12.4 9.9 12.4 9.9 13.4 13.4  9.9   -12.4  -9.9 -12.4]/10;
-                                 % diaphragma-array distance in [cm] (poloidal)
-   da2=[37 34.4 37 55.9 59.4 59.4 55.9 37 34.4 37]/10;
-                                 % dist to diaphragm in toroidal direction [cm];
-   deltada=[   -0.0311 0  -0.0458   -0.1179   -0.0615   -0.1105 ...
-			  -0.0510   -0.0515  0  -0.3223]/10;
-  deltada(4)=0;
-  deltada(6)=0;
-                                                                 
-
-   da=da+ deltada;
-
-   da2=da2+deltada;
-
-
-   d1=0.90;			 % detector width      in mm
-   d2=4.0;			 % detector length     in mm
-   b1=0.800;			 % aperture width      in mm (pol.)
-   b2=[8 8 8 15 15 15 15 8 8 8];			 
-				 % aperture length     in mm (tor.)
-   b3x=0.020;			 % aperture thickness  in mm (poloidal)
-   b3y=0;			 % aperture thickness  in mm (toroidal)
-end
-
-
-
-%======== calculation of the chords of view ===================================
-
-nact=sum(fans);
-iact=find(fans);
-ndet=20;
-ncam=10;
-
-
-
-% ---- apertures: ------------------
-
-xap=ones(ndet,1)*xpos(iact);
-xap=xap(:)';
-yap=ones(ndet,1)*ypos(iact);
-yap=yap(:)';
-
-% ---- detectors: ------------------
-
-vorz(find(vangle==90))=(-1)^(cw+1)*ones(size(find(vangle==90)));
-vorz(find(vangle==0))=(-1)^cw*ones(size(find(vangle==0)));
-vorz(find(vangle==-90))=(-1)^cw*ones(size(find(vangle==-90)));
-
-
-dete=(-9.025:0.950:9.025)'/10*vorz(iact)+ones(ndet,1)*ae(iact);
-
-dum_ae=dete(:)';
-
-dum_vangle=ones(ndet,1)*vangle(iact);
-dum_vangle=dum_vangle(:)';
-
-
-
-ivert=find(dum_vangle==90 | dum_vangle==-90);
-ihori=find(dum_vangle==0);
-
-dum_da=ones(ndet,1)*da(iact);
-dum_da=dum_da(:)';
-
-dxd=zeros(1,ndet*nact);
-dyd=zeros(1,ndet*nact);
-
-dxd(ivert)=dum_ae(ivert);
-dxd(ihori)=dum_da(ihori);
-
-dyd(ivert)=dum_da(ivert);
-dyd(ihori)=dum_ae(ihori);
-
-xdet=xap+dxd;
-ydet=yap+dyd;
-
-
-%plot_vessel(rzvin,rzvout)
-%hold on
-% plot(xap,yap,'.g',xdet,ydet,'.m')
-
-
-% ---- calculate the equations of lines of sight
-
-m=(ydet-yap)./(xdet-xap);
-b=ydet-m.*xdet;
-
-nl=length(xdet);
-xchord=zeros(2,nl);
-ychord=zeros(2,nl);
-
-
-xchord(1,:)=xdet;ychord(1,:)=ydet;
-
-iup=find(dum_vangle==90);
-isi=find(dum_vangle==0);
-ido=find(dum_vangle==-90);
-
-if ~isempty(iup)
-ychord(2,iup)=ymin*ones(size(iup));
-xchord(2,iup)=(ychord(2,iup)-b(iup))./m(iup);
-end
-if ~isempty(ido)
-ychord(2,ido)=ymax*ones(size(ido));
-xchord(2,ido)=(ychord(2,ido)-b(ido))./m(ido);
-end
-if ~isempty(isi)
-xchord(2,isi)=xmin*ones(size(isi));
-ychord(2,isi)=m(isi).*xchord(2,isi)+b(isi);
-end
-
-ileft=find(xchord(2,:)<xmin);
-
-if ~isempty(ileft)
-xchord(2,ileft)=xmin*ones(size(ileft));
-ychord(2,ileft)=m(ileft).*xchord(2,ileft)+b(ileft);
-end
-irig=find(xchord(2,:)>xmax);
-if ~isempty(irig)
-xchord(2,irig)=xmax*ones(size(irig));
-ychord(2,irig)=m(irig).*xchord(2,irig)+b(irig);
-end
-
-
-
-%======== prepare output ======================================================
-
-vangle=vangle(iact);
-
-%======== calculation of angular correction factors, if necessary =============
-
-if i_detec==2 & exist('angular_fact_2.mat')==2
-
-		disp('loading angular_fact_2')
-		load angular_fact_2
-
-
-elseif i_detec==1 & exist('angular_fact_1.mat')==2
-
-		disp('loading angular_fact_1')
-		load angular_fact_1
-
-else
-
-		aomega=zeros(ndet,ncam);
-		angfact=ones(ndet,ncam);
-
-
-		for l=1:sum(fans)
-
-%		Z0X=abs(da(iact(l))*10)
-%	 	Z0Y=abs(da2(iact(l))*10)
-%		X0=ae(iact(l))*10 % back to mm, sorry about that...
-%		X0=X0*vorz(iact(l))
-%		B2=b2(iact(l))
-%		AOMEGA=etendue_n(b1,B2,b3x,b3y,Z0X,Z0Y,X0,cw);
-
-		b1x=0.8;
-		b1y=6;
-		b1z=0.02;
-		b2x=10000000;
-		b2y=b2(iact(l));
-		b2z=0;
-		z01=abs(da(iact(l))*10);
-		z02=abs(da2(iact(l))*10);
-		X0=ae(iact(l))*10*vorz(iact(l));
-keyboard
-		AOMEGA=etendue_n2(b1x,b1y,b1z,b2x,b2y,b2z,z01,z02,X0,cw);
-
-			aomega(:,iact(l))=AOMEGA(:,1);
-
-		
-		end
-
-
-		indm=min(find(aomega==max(aomega(:))));
-		aomegan=aomega/aomega(indm);
-		nonz=find(aomega);
-		angfact(nonz)=ones(size(nonz))./aomegan(nonz);
-
-		angfact=round(1000*angfact)/1000;
-
-
-		unitstring='units aomega: mm^2 * sterad';
-
-		if i_detec==1
-			save angular_fact_1 angfact aomega unitstring
-		elseif i_detec==2
-			save angular_fact_2 angfact aomega unitstring
-		end
-
-end
-return
-
-th=atan(diff(ychord)./diff(xchord));
-thp=th;
-neg=find(thp<0);
-thp(neg)=180+thp(neg);
-thdet=ones(1,20*sum(fans));
-for k=1:sum(fans)
-     thdet((k-1)*20+1:k*20)=vangle(k)*ones(1,20);
-end
-angles=thdet-thp;
-mist=find(angles<0 & abs(angles)>90);
-angles(mist)=angles(mist)+180;
-th_inc=angles*pi/180;
-
-
-% ---- correct for the edges of tcv ( some chords may be too long )
-
-
-
-
-
-	down=find(xcont>xedge & ycont<-yedge);
-	up=find(xcont>xedge & ycont>yedge);
-	cd=polyfit(xcont(down),ycont(down),1);
-	cu=polyfit(xcont(up),ycont(up),1);
-
-
-	iu1=find(xchord(1,:)>xedge & ychord(1,:)>0 & dum_vangle==-90 );
-	if ~isempty(iu1)
-		xchord(1,iu1)=-(b(iu1)-cu(2))./(m(iu1)-cu(1)+eps);
-		ychord(1,iu1)=m(iu1).*xchord(1,iu1)+b(iu1);
-	end
-
-	iu2=find(xchord(2,:)>xedge & ychord(2,:)>0 & ychord(1,:) & ....
-			dum_vangle==-90);
-	if ~isempty(iu2)
-		xchord(2,iu2)=-(b(iu2)-cu(2))./(m(iu2)-cu(1)+eps);
-		ychord(2,iu2)=m(iu2).*xchord(2,iu2)+b(iu2);
-	end
-
-	id1=find(xchord(1,:)>xedge & ychord(1,:)<0 & dum_vangle==90);
-	if ~isempty(id1)
-		xchord(1,id1)=-(b(id1)-cd(2))./(m(id1)-cd(1)+eps);
-		ychord(1,id1)=m(id1).*xchord(1,id1)+b(id1);
-	end
-
-	id2=find(xchord(2,:)>xedge & ychord(2,:)<0 & dum_vangle==90);
-	if ~isempty(id2)
-		xchord(2,id2)=-(b(id2)-cd(2))./(m(id2)-cd(1)+eps);
-		ychord(2,id2)=m(id2).*xchord(2,id2)+b(id2);
-	end
-
-	ilow=find(ychord(1,:)<ymin);
-	ihig=find(ychord(1,:)>ymax);
-	ilef=find(xchord(1,:)<xmin);
-	irig=find(xchord(1,:)>xmax);
-	if ~isempty(ilow)
-		ychord(1,ilow)=ymin*ones(size(ilow));
-		xchord(1,ilow)=ymin./m(ilow)-b(ilow)./m(ilow);
-	end
-	if ~isempty(ihig)
-		ychord(1,ihig)=ymax*ones(size(ihig));
-		xchord(1,ihig)=ymax./m(ihig)-b(ihig)./m(ihig);
-	end
-	if ~isempty(ilef)
-		xchord(1,ilef)=xmin*ones(size(ilef));
-		ychord(1,ilef)=m(ilef)*xmin+b(ilef);
-	end
-	if ~isempty(irig)
-		xchord(1,irig)=xmax*ones(size(irig));
-		ychord(1,irig)=m(irig)*xmax+b(irig);
-	end
-
-
-	ilow=find(ychord(2,:)<ymin);
-	ihig=find(ychord(2,:)>ymax);
-	ilef=find(xchord(2,:)<xmin);
-	irig=find(xchord(2,:)>xmax);
-	if ~isempty(ilow)
-		ychord(2,ilow)=ymin*ones(size(ilow));
-		xchord(2,ilow)=ymin./m(ilow)-b(ilow)./m(ilow);
-	end
-	if ~isempty(ihig)
-		ychord(2,ihig)=ymax*ones(size(ihig));
-		xchord(2,ihig)=ymax./m(ihig)-b(ihig)./m(ihig);
-	end
-	if ~isempty(ilef)
-		xchord(2,ilef)=xmin*ones(size(ilef));
-		ychord(2,ilef)=m(ilef)*xmin+b(ilef);
-	end
-	if ~isempty(irig)
-		xchord(2,irig)=xmax*ones(size(irig));
-		ychord(2,irig)=m(irig)*xmax+b(irig);
-	end
-
-
-
-
-