diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 41907f55ceb05140887dd8ef60ba518361cf8774..4680896b3ffed34a524be28c15d1f0e7d3cd711f 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -56,7 +56,7 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(shot,data_req
 %    a5=gdat(48836,'ip'); % standard call
 %    a6=gdat(48836,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (note all lowercase in output)
 %    a7 = gdat(48836,'static("r_m")[$1]'',''001'); % note first and last quote of tdi argument added by gdat
-%    a7 = gdat(48836,'tcv_eq(''''psi'''',''''liuqe.m'''')'); % do not use double quote char so 'liuqe',11 will work
+%    a7 = gdat(48836,'tcv_eq("psi","liuqe.m")'); % do not use double quote char so 'liuqe',11 will work
 %    a8 = gdat(64770,['\magnetics::bpol_003[*,$1]'',{''001''; ''030''}']); one string full expression
 %    a9 = gdat(64770,{'\magnetics::vloop[*,$1]',{'''001''','''A_001'''}}); cell array (other option)
 %
@@ -230,6 +230,7 @@ end
 
 % extract parameters from pairs of varargin:
 i_liuqe_set_in_pairs = 0; % to know if liuqe was not set explicitely thus use value in expression later
+if isstruct(data_request) && isfield(data_request,'liuqe'), i_liuqe_set_in_pairs = 1; end % given in input parameter structure
 if (nargin>=ivarargin_first_char)
   if mod(nargin-ivarargin_first_char+1,2)==0
     for i=1:2:nargin-ivarargin_first_char+1
@@ -290,12 +291,35 @@ shot = gdat_data.shot;
 data_request_eff = gdat_data.gdat_params.data_request;
 error_status = 6; % at least reached this level
 
+mapping_for_tcv = tcv_requests_mapping(data_request_eff,shot);
+gdat_data.label = mapping_for_tcv.label;
+
 liuqe_version = 1;
-if isfield(gdat_data.gdat_params,'liuqe') && ~isempty(gdat_data.gdat_params.liuqe)
+if isfield(gdat_data.gdat_params,'liuqe') && ~isempty(gdat_data.gdat_params.liuqe) && i_liuqe_set_in_pairs==1
   liuqe_version = gdat_data.gdat_params.liuqe;
 else
-  gdat_data.gdat_params.liuqe = liuqe_version;
+  % if no specific liuqe requested in the parameter option, but specified e.g. in tcv_eq, use that one otherwise default
+  if ~iscell(mapping_for_tcv.expression)
+    bb{1} = mapping_for_tcv.expression;
+  else
+    bb = mapping_for_tcv.expression;
+  end
+  for i=1:numel(bb)
+    if any(regexpi(bb{i},'fbte','once'))
+      liuqe_version = -1
+    else
+      ij = regexpi(bb{i},'LIUQE\.M.?','once');
+      if ~isempty(ij) && any(regexpi(bb{i},'LIUQE\.M[2,3]','once'))
+        liuqe_version = str2num(bb{i}(ij+7))
+      end
+      ij = regexpi(bb{i},'LIUQE[^\.]','once');
+      if ~isempty(ij) && any(regexpi(bb{i},'LIUQE[2,3]','once'))
+        liuqe_version = str2num(bb{i}(ij+5))
+      end
+    end
+  end
 end
+gdat_data.gdat_params.liuqe = liuqe_version;
 liuqe_matlab = 1; % now default should be matlab liuqe nodes
 if liuqe_version<0 || (liuqe_version > 10 && liuqe_version < 20)
   liuqe_matlab = 0;
@@ -310,8 +334,6 @@ if liuqe_version_eff==2 || liuqe_version_eff==3
   substr_liuqe_tcv_eq = num2str(liuqe_version_eff);
 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 >=900'000
 begstr = '';
 if (iscell(mapping_for_tcv.expression) || isempty(strfind(mapping_for_tcv.expression,'\rtc::'))) && ...
@@ -332,6 +354,8 @@ switch liuqe_version
  case {12,13}, liuqe_ext=['_' num2str(mod(liuqe_version,10))]; psitbx_str=['LIUQE' num2str(mod(liuqe_version,10))];liuqe_matlab = 0;
  otherwise, error(['Unknown LIUQE version, liuqe = ' num2str(liuqe_version)]);
 end
+% entry in gdat_request superseded by 'liuqe' option (so remove replacement later)
+mapping_for_tcv.expression = regexprep(mapping_for_tcv.expression,'FBTE([",''])|LIUQE[2-3]?([",''])|LIUQE.M[2-3]?([",''])',[psitbx_str '$1'],'preservecase');
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -411,22 +435,14 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi')
       return
     end
   else
-    do_liuqem2liuqef = 0;
+    do_liuqem2liuqef = 1;
     if liuqe_version_eff==-1
       % trialindx not relevant for fbte
       mapping_for_tcv_expression_eff = mapping_for_tcv.expression;
       if length(mapping_for_tcv.expression)>8 && strcmp(lower(mapping_for_tcv.expression(1:8)),'\results')
         mapping_for_tcv_expression_eff = mapping_for_tcv.expression(11:end);
-      elseif findstr('tcv_eq',lower(mapping_for_tcv.expression))
-        ij = regexpi(mapping_for_tcv.expression,'''''');
-        if length(ij)<2
-          disp(['expected more double quotes in mapping_for_tcv.expression = ' mapping_for_tcv.expression]);
-          return
-        else
-          mapping_for_tcv_expression_eff = mapping_for_tcv.expression(ij(1)+2:ij(2)-1);
-        end
       end
-      eval_expr = ['tdi(''' begstr mapping_for_tcv_expression_eff substr_liuqe ''');']
+      eval_expr = ['tdi(''' begstr mapping_for_tcv_expression_eff substr_liuqe ''');'];
     else
       aaa_tmp = regexp(mapping_for_tcv.expression,{'toray','conf','ibs','proffit','astra'});
       if ~isempty(gdat_data.gdat_params.trialindx) && gdat_data.gdat_params.trialindx >= 0 ...
@@ -446,66 +462,25 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi')
         if ~isempty(ij)
           mapping_for_tcv.expression(ij+5:ij+6) = substr_liuqe;
         end
-        ij = regexpi(mapping_for_tcv.expression,'LIUQE\.M.?','once');
-        if ~isempty(ij)
-          ichar_after_liuqe = 7;
-          if strcmp(mapping_for_tcv.expression(ij+ichar_after_liuqe),'2') || ...
-                strcmp(mapping_for_tcv.expression(ij+ichar_after_liuqe),'3')
-            if i_liuqe_set_in_pairs==0
-              gdat_params.liuqe = str2num(mapping_for_tcv.expression(ij+ichar_after_liuqe));
-              gdat_data.gdat_params = gdat_params;
-              substr_liuqe_tcv_eq = mapping_for_tcv.expression(ij+ichar_after_liuqe);
-            end
-            ichar_after_liuqe = 8;
-          end
-          mapping_for_tcv.expression = [mapping_for_tcv.expression(1:ij+6) substr_liuqe_tcv_eq mapping_for_tcv.expression(ij+ichar_after_liuqe:end)];
-        else
-          % check if liuqe, liuqe2 or liuqe3 is given in expression
-          ij = regexpi(mapping_for_tcv.expression,'LIUQE[^\.]','once');
-          if ~isempty(ij)
-            ichar_after_liuqe = 5;
-            if strcmp(mapping_for_tcv.expression(ij+ichar_after_liuqe),'2') || ...
-                  strcmp(mapping_for_tcv.expression(ij+ichar_after_liuqe),'3')
-              if i_liuqe_set_in_pairs==0
-                gdat_params.liuqe = 10+str2num(mapping_for_tcv.expression(ij+ichar_after_liuqe));
-                gdat_data.gdat_params = gdat_params;
-                substr_liuqe_tcv_eq = mapping_for_tcv.expression(ij+ichar_after_liuqe);
-              end
-              ichar_after_liuqe = 6;
-            else
-              if i_liuqe_set_in_pairs==0
-                gdat_params.liuqe = 11;
-                gdat_data.gdat_params = gdat_params;
-                substr_liuqe_tcv_eq = '';
-              end
-            end
-            if i_liuqe_set_in_pairs==1 && liuqe_matlab==1
-              % enforce matlab liuqe version even matlab if asked for
-              substr_liuqe_tcv_eq = ['.M' substr_liuqe_tcv_eq];
-            end
-            mapping_for_tcv.expression = [mapping_for_tcv.expression(1:ij+4) substr_liuqe_tcv_eq mapping_for_tcv.expression(ij+ichar_after_liuqe:end)];
-          end
-        end
-      else
-        ij = regexpi(mapping_for_tcv.expression,'LIUQE\.M.?','once');
-        if ~isempty(ij)
-          mapping_for_tcv.expression = [mapping_for_tcv.expression(1:ij+4) substr_liuqe_tcv_eq mapping_for_tcv.expression(ij+7:end)];
-          do_liuqem2liuqef = 1;
-        end
       end
       eval_expr = ['tdi(''' mapping_for_tcv.expression ''');'];
     end
     % special cases for matching liuqe fortran and matlab
     if liuqe_matlab == 0 && do_liuqem2liuqef==1
-      % assume tcv_eq(''keyword'',''LIUQE..'')
+      % assume tcv_eq(''keyword'',''LIUQE..'') or tcv_eq("keyword","LIUQE..")
       liuqe_fortran_matlab = liuqefortran2liuqematlab;
-      ij = regexpi(eval_expr,'''''');
-      if numel(ij)>=2
-        liuqe_keyword = eval_expr(ij(1)+2:ij(2)-1);
+      eval_expr_eff = strrep(eval_expr,'''','');
+      eval_expr_eff = strrep(eval_expr_eff,'"','');
+      ij = regexpi(eval_expr_eff,',');
+      if isempty(ij)
+        disp(['expected tdi(tcv_eq(...,...) with quotes or double quotes in eval_expr = ' eval_expr]);
+        return
+      else
+        liuqe_keyword = eval_expr_eff(4+8:ij(1)-1);
       end
       ij_row=strmatch(liuqe_keyword,liuqe_fortran_matlab(:,2),'exact');
       if ~isempty(ij_row)
-        eval_expr = [eval_expr(1:ij(1)+1) liuqe_fortran_matlab{ij_row,1} eval_expr(ij(1)+2+length(liuqe_keyword):end)];
+        eval_expr = regexprep(eval_expr,['\<' liuqe_keyword '\>'],liuqe_fortran_matlab{ij_row,1},'preservecase');
       end
     end
     aatmp=eval(eval_expr);
@@ -525,7 +500,7 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi')
       end
     end
   end
-  if isempty(aatmp.data) || (isempty(aatmp.dim) && ischar(aatmp.data) && ~isempty(strfind(lower(aatmp.data),'no data')))% || ischar(aatmp.data) (to add?)
+  if isempty(aatmp.data) || (isempty(aatmp.dim) && ischar(aatmp.data) && ~isempty(strfind(lower(aatmp.data),'no data')) ) % || ischar(aatmp.data) (to add?)
     if (gdat_params.nverbose>=1); warning(['problems loading data for ' eval_expr ' for data_request= ' data_request_eff]); end
     if (gdat_params.nverbose>=3); disp('check .gdat_request list'); end
     return
@@ -688,7 +663,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     if strcmp(data_request_eff,'r_geom_rho'); data_request_eff = 'rgeom_rho'; end
     % compute average minor or major radius (on z=zaxis normally)
     if liuqe_matlab==0
-      nodenameeff=['tcv_eq(''r_max_psi'',''LIUQE' substr_liuqe_tcv_eq ''')'];
+      nodenameeff=['tcv_eq("r_max_psi","' psitbx_str '")'];
       rmaxpsi=tdi(nodenameeff);
       ijnan = find(isnan(rmaxpsi.data));
       if isempty(rmaxpsi.data) || isempty(rmaxpsi.dim) || ischar(rmaxpsi.data) || ...
@@ -697,7 +672,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
         return
       end
-      nodenameeff2=['tcv_eq(''r_min_psi'',''LIUQE' substr_liuqe_tcv_eq ''')'];
+      nodenameeff2=['tcv_eq("r_min_psi","' psitbx_str '")'];
       rminpsi=tdi(nodenameeff2);
       ijnan = find(isnan(rminpsi.data));
       if isempty(rminpsi.data) || isempty(rminpsi.dim) || ...
@@ -730,7 +705,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         gdat_data.dim = rmaxpsi.dim;
         gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
         gdat_data.x = gdat_data.dim{1};
-        gdat_data.dimunits = rmaxpsi.dimunits(2);
+        gdat_data.dimunits = rmaxpsi.dimunits;
       else
         gdat_data.dim = rmaxpsi.dim(2);
         gdat_data.t = gdat_data.dim{mapping_for_tcv.timedim};
@@ -742,19 +717,19 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     else
       if isempty(substr_liuqe); substr_liuqe = '_1'; end
       if strcmp(data_request_eff,'a_minor') % to update when ready, still buggy, this should be rho,t
-        nodenameeff=['tcv_eq(''a_minor_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
+        nodenameeff=['tcv_eq("a_minor_mid","' psitbx_str '")'];
       elseif strcmp(data_request_eff,'a_minor_rho')
-        nodenameeff=['tcv_eq(''r_out_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
-        nodenameeff2=['tcv_eq(''r_in_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
-        nodenameeff=['tcv_eq(''r_out'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
-        nodenameeff2=['tcv_eq(''r_in'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
+        nodenameeff=['tcv_eq("r_out_mid","' psitbx_str '")'];
+        nodenameeff2=['tcv_eq("r_in_mid","' psitbx_str '")'];
+        nodenameeff=['tcv_eq("r_out","' psitbx_str '")'];
+        nodenameeff2=['tcv_eq("r_in","' psitbx_str '")'];
       elseif strcmp(data_request_eff,'rgeom')
-        nodenameeff=['tcv_eq(''r_geom_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
+        nodenameeff=['tcv_eq("r_geom_mid","' psitbx_str '")'];
       elseif strcmp(data_request_eff,'rgeom_rho')
-        nodenameeff=['tcv_eq(''r_out_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
-        nodenameeff2=['tcv_eq(''r_in_mid'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
-        nodenameeff=['tcv_eq(''r_out'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
-        nodenameeff2=['tcv_eq(''r_in'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
+        nodenameeff=['tcv_eq("r_out_mid","' psitbx_str '")'];
+        nodenameeff2=['tcv_eq("r_in_mid","' psitbx_str '")'];
+        nodenameeff=['tcv_eq("r_out","' psitbx_str '")'];
+        nodenameeff2=['tcv_eq("r_in","' psitbx_str '")'];
       else
         if (gdat_params.nverbose>=1); warning(['should not be in this case with data_request_eff = ' data_request_eff]); end
         return
@@ -795,10 +770,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.gdat_params.gdat_request = gdat_data.gdat_request;
     gdat_params = gdat_data.gdat_params;
     if liuqe_matlab==0
-      nodenameeff=['tcv_eq(''z_contour'',''LIUQE' substr_liuqe_tcv_eq ''')'];
+      nodenameeff=['tcv_eq("z_contour","' psitbx_str '")'];
     else
       if isempty(substr_liuqe); substr_liuqe = '_1'; end
-      nodenameeff=['tcv_eq(''z_edge'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
+      nodenameeff=['tcv_eq("z_edge","' psitbx_str '")'];
     end
     zcontour=tdi(nodenameeff);
     if isempty(zcontour.data) || isempty(zcontour.dim)  % || ischar(zcontour.data) (to add?)
@@ -850,10 +825,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       else
         if liuqe_matlab==0
           added_correction = 1.0; % correction already in liuqe.f
-          nodenameeff = ['tcv_eq(''BZERO'',''LIUQE' substr_liuqe_tcv_eq ''')'];
+          nodenameeff = ['tcv_eq("BZERO","' psitbx_str '")'];
         else
           if isempty(substr_liuqe); substr_liuqe = '_1'; end
-          nodenameeff=['tcv_eq(''BZERO'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
+          nodenameeff=['tcv_eq("BZERO","' psitbx_str '")'];
           added_correction = 0.996; % correction removed in liuqe.m at this stage
           added_correction_str = [num2str(added_correction) ' * '];
         end
@@ -1307,7 +1282,11 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.gas_request_flux = gasrequest;
     gdat_data.gas_request_flux.data = max(0.,72.41.*(gasrequest.data-0.6879)./(4.2673-0.6879));
     gdat_data.gas_request_flux.units = gasflux.units;
-    params_eff.data_request = '\draw_feedfor_gas:alim_001';
+    if shot < 78667
+      params_eff.data_request = '\draw_feedfor_gas:alim_001';
+    else
+      params_eff.data_request = '\draw_refs_gas:ref_001';
+    end
     gasrequest_ref = gdat_tcv(gdat_data.shot,params_eff); gasrequest_ref.units = 'V';
     gdat_data.gas_feedforward_volt = gasrequest_ref;
     gdat_data.gas_feedforward_flux = gasrequest_ref;
@@ -2727,9 +2706,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
    case {'q_rho'}
     % q profile on psi from liuqe
     if liuqe_matlab==0
-      nodenameeff=['tcv_eq(''q_psi'',''LIUQE' substr_liuqe_tcv_eq ''')'];
+      nodenameeff=['tcv_eq("q_psi","' psitbx_str '")'];
     else
-      nodenameeff=['tcv_eq(''q'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
+      nodenameeff=['tcv_eq("q","' psitbx_str '")'];
     end
     if liuqe_version_eff==-1
       nodenameeff=[begstr 'q_psi' substr_liuqe];
@@ -2766,10 +2745,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
    case {'rbphi_rho', 'rbtor_rho'}
     % R*Bphi(rho,t) from F from FFprime
     if liuqe_matlab==0
-      disp('not yet implemented for liuqe fortran')
+      fprintf('\n***** not yet implemented for liuqe fortran ******\n');
       return
     else
-      nodenameeff=['tcv_eq(''rbtor_rho'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
+      nodenameeff=['tcv_eq("rbtor_rho","' psitbx_str '")'];
     end
     if liuqe_version_eff==-1
       disp('not yet implemented for fbte')
@@ -2811,7 +2790,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     % need to get q_rho but to avoid loop for rhotor in grids_1d, get q_rho explicitely here
     params_eff = gdat_data.gdat_params;
     if liuqe_matlab==1
-      nodenameeff = ['tcv_eq(''tor_flux_tot'',''' psitbx_str ''')'];
+      nodenameeff = ['tcv_eq("tor_flux_tot","' psitbx_str '")'];
       phi_tor = tdi(nodenameeff);
     else
       phi_tor.data = [];
@@ -2819,7 +2798,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
     if ~any(any(isfinite(phi_tor.data))) || ischar(phi_tor.data)
       % no phi_tor, compute it from q profile
-      q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']);
+      q_rho = tdi(['tcv_eq("' liuqefortran2liuqematlab('q',1,liuqe_matlab) '","' psitbx_str '")']);
       if liuqe_matlab==0
         q_rho.dim{1} = sqrt(linspace(0,1,numel(q_rho.dim{1})));
       end
@@ -2847,7 +2826,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       if any(any(~isfinite(phi_tor.data(1:end-1,:))))
         only_edge = 1;
       end
-      q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']);
+      q_rho = tdi(['tcv_eq("' liuqefortran2liuqematlab('q',1,liuqe_matlab) '","' psitbx_str '")']);
       if liuqe_matlab==0
         q_rho.dim{1} = sqrt(linspace(0,1,numel(q_rho.dim{1})));
       end
@@ -2910,14 +2889,14 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'pprime', 'pprime_rho'}
     if liuqe_matlab==0
-      nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('pprime_rho',1,0) '",''' psitbx_str ''')'];
+      nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('pprime_rho',1,0) '","' psitbx_str '")'];
       tracetdi=tdi(nodenameeff);
       if ~isempty(tracetdi.dim) && length(tracetdi.dim)>=1
         tracetdi.dim{1} = sqrt(tracetdi.dim{1}); % correct x-axis psi_norm to rhopol
       end
       tracetdi.data = tracetdi.data ./2 ./pi; % correct node assumption (same for liuqe_fortran and fbte)
     else
-      nodenameeff = ['tcv_eq(''pprime_rho'',''' psitbx_str ''')'];
+      nodenameeff = ['tcv_eq("pprime_rho","' psitbx_str '")'];
       tracetdi=tdi(nodenameeff);
     end
     gdat_data.data = tracetdi.data;
@@ -2937,14 +2916,14 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
    case {'pressure', 'pressure_rho', 'p_rho'}
     if liuqe_matlab==0
-      nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('p_rho',1,0) '",''' psitbx_str ''')'];
+      nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('p_rho',1,0) '","' psitbx_str '")'];
       tracetdi=tdi(nodenameeff);
       if ~isempty(tracetdi.dim) && length(tracetdi.dim)>=1
         tracetdi.dim{1} = sqrt(tracetdi.dim{1}); % correct x-axis psi_norm to rhopol
       end
       tracetdi.data = tracetdi.data ./2 ./pi; % correct node assumption (same for liuqe_fortran and fbte)
     else
-      nodenameeff = ['tcv_eq(''p_rho'',''' psitbx_str ''')'];
+      nodenameeff = ['tcv_eq("p_rho","' psitbx_str '")'];
       tracetdi=tdi(nodenameeff);
     end
     gdat_data.data = tracetdi.data;
@@ -3008,7 +2987,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       tracetdi=tdi(nodenameeff);
       gdat_data.request_description = 'NaNs when smaller nb of boundary points at given time, can use \results::npts_contour';
     else
-      nodenameeff = ['tcv_eq(''' data_request_eff(1) '_edge'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
+      nodenameeff = ['tcv_eq("' data_request_eff(1) '_edge","' psitbx_str '")'];
       tracetdi=tdi(nodenameeff);
       gdat_data.request_description = 'lcfs on same nb of points for all times';
     end
@@ -3033,7 +3012,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     % need to get q_rho but to avoid loop for rhotor in grids_1d, get q_rho explicitely here
     params_eff = gdat_data.gdat_params;
     if liuqe_matlab==0
-      nodenameeff=['tcv_eq(''q_psi'',''LIUQE' substr_liuqe_tcv_eq ''')'];
+      nodenameeff=['tcv_eq("q_psi","' psitbx_str '")'];
       if liuqe_version_eff==-1
         nodenameeff=[begstr 'q_psi' substr_liuqe];
       end
@@ -3176,7 +3155,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         return
       end
     else
-      nodenameeff=['tcv_eq(''vol'',''' psitbx_str ''')'];
+      nodenameeff=['tcv_eq("vol","' psitbx_str '")'];
     end
     if liuqe_version_eff==-1
       data_request_eff = 'volume'; % only LCFS
@@ -3464,14 +3443,14 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ttprime', 'ttprime_rho'}
     if liuqe_matlab==0
-      nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('ttprime_rho',1,0) '",''' psitbx_str ''')'];
+      nodenameeff = ['tcv_eq("' liuqefortran2liuqematlab('ttprime_rho',1,0) '","' psitbx_str '")'];
       tracetdi=tdi(nodenameeff);
       if ~isempty(tracetdi.dim) && length(tracetdi.dim)>=1
         tracetdi.dim{1} = sqrt(tracetdi.dim{1}); % correct x-axis psi_norm to rhopol
       end
       tracetdi.data = tracetdi.data ./2 ./pi; % correct node assumption (same for liuqe_fortran and fbte)
     else
-      nodenameeff = ['tcv_eq(''ttprime_rho'',''' psitbx_str ''')'];
+      nodenameeff = ['tcv_eq("ttprime_rho","' psitbx_str '")'];
       tracetdi=tdi(nodenameeff);
     end
     gdat_data.data = tracetdi.data;
diff --git a/matlab/TCV/tcv_requests_mapping.m b/matlab/TCV/tcv_requests_mapping.m
index 47be410cad59e6f6fdceca56a4e9ae9c05e55821..1bcfcdd585d6c136231618b0c2cdf8598272e301 100644
--- a/matlab/TCV/tcv_requests_mapping.m
+++ b/matlab/TCV/tcv_requests_mapping.m
@@ -53,17 +53,17 @@ switch lower(data_request)
   mapping.timedim = 1;
   mapping.label = 'area';
   mapping.method = 'tdiliuqe';
-  mapping.expression = 'tcv_eq(''''area'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("area","LIUQE.M")';
  case 'area_rho'
   mapping.timedim = 2;
   mapping.label = 'area_rho';
   mapping.method = 'tdiliuqe';
-  mapping.expression = 'tcv_eq(''''area_rho'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("area_rho","LIUQE.M")';
  case 'area_edge'
   mapping.timedim = 1;
   mapping.label = 'area\_lcfs';
   mapping.method = 'tdiliuqe';
-  mapping.expression = 'tcv_eq(''''area_edge'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("area_edge","LIUQE.M")';
  case 'b0'
   mapping.timedim = 1;
   mapping.label = 'B_0';
@@ -74,7 +74,7 @@ switch lower(data_request)
   mapping.label = '\beta';
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::beta_tor';
-  mapping.expression = 'tcv_eq(''''beta_tor'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("beta_tor","LIUQE.M")';
  case {'betan', 'beta_tor_norm'}
   mapping.timedim = 1;
   mapping.label = '\beta_N';
@@ -86,7 +86,7 @@ switch lower(data_request)
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::beta_pol';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:beta_pol';
-  mapping.expression = 'tcv_eq(''''beta_pol'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("beta_pol","LIUQE.M")';
  case 'cxrs'
   mapping.timedim = 2;
   mapping.label = 'cxrs';
@@ -102,7 +102,7 @@ switch lower(data_request)
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::delta_edge';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:delta_edge';
-  mapping.expression = 'tcv_eq(''''delta_edge'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("delta_edge","LIUQE.M")';
   % mapping.method = 'expression';
   % mapping.expression = ['tdi(''\results::q_psi'');']; % for tests
  case 'delta_bottom'
@@ -111,19 +111,19 @@ switch lower(data_request)
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::delta_ed_bot';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:delta_bot';
-  mapping.expression = 'tcv_eq(''''delta_ed_bot'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("delta_ed_bot","LIUQE.M")';
  case 'delta_rho'
   mapping.timedim = 2;
   mapping.method = 'tdiliuqe';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:delta';
-  mapping.expression = 'tcv_eq(''''delta'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("delta","LIUQE.M")';
  case 'delta_top'
   mapping.timedim = 1;
   mapping.label = 'delta\_top';
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::delta_ed_top';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:delta_top';
-  mapping.expression = 'tcv_eq(''''delta_ed_top'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("delta_ed_top","LIUQE.M")';
  case 'ece'
   mapping.timedim = 2;
   mapping.method = 'switchcase';
@@ -204,24 +204,40 @@ switch lower(data_request)
   mapping.label = 'Plasma current';
   mapping.timedim = 1;
   mapping.method = 'tdiliuqe';
-  mapping.expression = 'tcv_eq(''''i_pl'''',''''LIUQE.M'''')'; % to be able to get ip consistent with relevant LIUQE value
+  mapping.expression = 'tcv_eq("i_pl","LIUQE.M")'; % to be able to get ip consistent with relevant LIUQE value
+ case 'i_pol'
+  disp('use ''ipol'' to get values from magnetics');
+  mapping.timedim = 2;
+  mapping.method = 'expression';
+  mapping.label = 'ipol from liuqe';
+  mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''tcv_eq("i_pol","LIUQE.M")''; ' ...
+                    'gdat_tmp=gdat_tcv([],params_eff);coilsname=mdsvalue(''dim_of(\magnetics::ipol,1)''); ' ...
+                    'gdat_tmp.dimunits{1}=coilsname([3:18 20 1:2]);'];
+ case 'ipol'
+  disp('use ''i_pol'' to get values from liuqe');
+  mapping.timedim = 2; % changed to match usual time as last one and match liuqe i_pol
+  mapping.method = 'expression';
+  mapping.label = 'ipol from magnetics';
+  mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''\magnetics::ipol''; ' ...
+                    'gdat_tmp=gdat_tcv([],params_eff);gdat_tmp.dimunits{2}=''s'';gdat_tmp.dimunits{1}=gdat_tmp.dim{2};' ...
+                    'gdat_tmp.dim{2}=gdat_tmp.dim{1};gdat_tmp.dim{1}=[0:numel(gdat_tmp.dimunits{1})-1]'';gdat_tmp.data=gdat_tmp.data'''];
  case 'kappa'
   mapping.timedim = 1;
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::kappa_edge';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:kappa_edge';
-  mapping.expression = 'tcv_eq(''''kappa_edge'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("kappa_edge","LIUQE.M")';
  case 'kappa_rho'
   mapping.timedim = 2;
   mapping.method = 'tdiliuqe';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:kappa';
-  mapping.expression = 'tcv_eq(''''kappa'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("kappa","LIUQE.M")';
  case 'li'
   mapping.timedim = 1;
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::l_i';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:l_i_3';
-  mapping.expression = 'tcv_eq(''''l_i_3'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("l_i_3","LIUQE.M")';
  case 'mhd'
   mapping.timedim = 1;
   mapping.label = 'n=1,2, etc';
@@ -287,7 +303,7 @@ switch lower(data_request)
  case {'phi_tor', 'phitor', 'toroidal_flux'}
   mapping.timedim = 2;
   mapping.method = 'tdiliuqe';
-  mapping.expression = 'tcv_eq(''''tor_flux_tot'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("tor_flux_tot","LIUQE.M")';
   % node not filled in and in any case need special case for
   mapping.label = 'toroidal_flux';
   mapping.method = 'switchcase';
@@ -318,7 +334,7 @@ switch lower(data_request)
   mapping.timedim = 1;
   mapping.label = 'psi(R,Z)';
   mapping.method = 'tdiliuqe';
-  mapping.expression = 'tcv_eq(''''psi'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("psi","LIUQE.M")';
  case 'powers'
   mapping.timedim = 1;
   mapping.label = 'various powers';
@@ -336,7 +352,7 @@ switch lower(data_request)
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::psi_axis';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:psi_axis';
-  mapping.expression = 'tcv_eq(''''psi_axis'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("psi_axis","LIUQE.M")';
   mapping.label = 'psi\_axis with psi_edge=0';
  case 'psi_edge'
   mapping.timedim = 1;
@@ -345,25 +361,25 @@ switch lower(data_request)
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::surface_flux';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:psi_surf';
-  mapping.expression = 'tcv_eq(''''psi_surf'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("psi_surf","LIUQE.M")';
  case 'q0'
   mapping.timedim = 1;
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::q_zero';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:q_axis';
-  mapping.expression = 'tcv_eq(''''q_axis'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("q_axis","LIUQE.M")';
  case 'q95'
   mapping.timedim = 1;
   mapping.method = 'tdiliuqe';
   % mapping.expression = '\results::q_95';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:q_95';
-  mapping.expression = 'tcv_eq(''''q_95'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("q_95","LIUQE.M")';
  case 'qedge'
   mapping.timedim = 1;
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::q_edge';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:q_edge';
-  mapping.expression = 'tcv_eq(''''q_edge'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("q_edge","LIUQE.M")';
  case 'q_rho'
   mapping.timedim = 2;
   mapping.label = 'q';
@@ -377,19 +393,19 @@ switch lower(data_request)
 % $$$   mapping.method = 'tdiliuqe';
 % $$$   mapping.expression = '\results::r_contour';
 % $$$   mapping.expression = '\tcv_shot::top.results.equil_1.results:r_rho'; % several flux surfaces R coordinates (irho,itheta,t)
-% $$$   mapping.expression = 'tcv_eq(''''r_rho'''',''''LIUQE.M'''')';
+% $$$   mapping.expression = 'tcv_eq("r_rho","LIUQE.M")';
  case 'r_contour_edge'
   mapping.timedim = 2;
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::r_contour';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:r_surf'; % LCFS R coordinates (r,t)
-  mapping.expression = 'tcv_eq(''''r_edge'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("r_edge","LIUQE.M")';
 
   mapping.label = 'R\_lcfs';
   mapping.method = 'switchcase';
 
 % $$$   mapping.method = 'expression';
-% $$$   mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''tcv_eq(''''''''r_edge'''''''',''''''''LIUQE.M'''''''')''; ' ...
+% $$$   mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''tcv_eq("r_edge","LIUQE.M")''; ' ...
 % $$$                     'gdat_tmp=gdat_tcv(shot,params_eff);gdat_tmp.data=gdat_tmp.data(:,:,end);' ...
 % $$$                     'gdat_tmp.dim=gdat_tmp.dim(1:2);gdat_tmp.dimunits=gdat_tmp.dimunits(1:2);'];
  case {'rgeom', 'r_geom'}
@@ -422,7 +438,7 @@ switch lower(data_request)
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::r_axis';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:r_axis';
-  mapping.expression = 'tcv_eq(''''r_axis'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("r_axis","LIUQE.M")';
  case {'rtc','scd'}
   if strcmp('scd',lower(data_request)); error('********* should call with request****** ''rtc'''); end
   mapping.timedim = 1;
@@ -433,13 +449,13 @@ switch lower(data_request)
   mapping.timedim = 1;
   mapping.label = 'surface\_tor(rho)';
   mapping.method = 'tdiliuqe';
-  mapping.expression = 'tcv_eq(''''surf'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("surf","LIUQE.M")';
   mapping.help = 'toroidal surface'; % works for method tdi..., in method expression, gdat_tmp.help should be defined
  case {'surface', 'surface_edge'}
   mapping.timedim = 1;
   mapping.label = 'surface';
   mapping.method = 'expression';
-  mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''tcv_eq(''''''''surf'''''''',''''''''LIUQE.M'''''''')'';' ...
+  mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''tcv_eq("surf","LIUQE.M")'';' ...
                     'gdat_tmp=gdat_tcv(shot,params_eff); gdat_tmp.dim = {gdat_tmp.t}; gdat_tmp.x=[]; gdat_tmp.data= gdat_tmp.data(end,:);' ...
                     'gdat_tmp.dimunits{1}=''s'';gdat_tmp.help=''toroidal surface of LCFS'';'];
  case 'sxr'
@@ -508,7 +524,7 @@ switch lower(data_request)
 % $$$     corr = [4e3,4.5e3];
 % $$$     if shot==57732; corr=[4650,4650]; end
 % $$$     mapping.method = 'expression';
-% $$$     mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''tcv_eq(''''''''w_mhd'''''''',''''''''LIUQE.M'''''''')'';' ...
+% $$$     mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''tcv_eq("w_mhd","LIUQE.M")'';' ...
 % $$$                     'gdat_tmp=gdat_tcv(shot,params_eff);ij=find(gdat_tmp.t>0.5&gdat_tmp.t<1.03);' ...
 % $$$                     'aa=interp1([' num2str(time_for_corr(1)) ' ' num2str(time_for_corr(2)) ...
 % $$$                     '],[' num2str(corr(1)) ' ' num2str(corr(2)) '],gdat_tmp.t(ij));' ...
@@ -517,26 +533,26 @@ switch lower(data_request)
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::total_energy';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:w_mhd';
-  mapping.expression = 'tcv_eq(''''w_mhd'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("w_mhd","LIUQE.M")';
 % $$$   end
 % $$$  case 'z_contour'   % z_rho not yet implemented
 % $$$   mapping.timedim = 2;
 % $$$   mapping.method = 'tdiliuqe';
 % $$$   mapping.expression = '\results::z_contour';
 % $$$   mapping.expression = '\tcv_shot::top.results.equil_1.results:z_rho'; % several flux surfaces Z coordinates (irho,itheta,t)
-% $$$   mapping.expression = 'tcv_eq(''''z_rho'''',''''LIUQE.M'''')';
+% $$$   mapping.expression = 'tcv_eq("z_rho","LIUQE.M")';
  case 'z_contour_edge'
   mapping.timedim = 2;
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::z_contour';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:z_surf'; % LCFS Z coordinates (r,t)
-  mapping.expression = 'tcv_eq(''''z_edge'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("z_edge","LIUQE.M")';
 
   mapping.label = 'Z\_lcfs';
   mapping.method = 'switchcase';
 
 % $$$   mapping.method = 'expression';
-% $$$   mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''tcv_eq(''''''''z_edge'''''''',''''''''LIUQE.M'''''''')''; ' ...
+% $$$   mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''tcv_eq("z_edge","LIUQE.M")''; ' ...
 % $$$                     'gdat_tmp=gdat_tcv(shot,params_eff);gdat_tmp.data=gdat_tmp.data(:,:,end);' ...
 % $$$                     'gdat_tmp.dim=gdat_tmp.dim(1:2);gdat_tmp.dimunits=gdat_tmp.dimunits(1:2);'];
  case 'zeff'
@@ -554,7 +570,7 @@ switch lower(data_request)
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::z_axis';
   mapping.expression = '\tcv_shot::top.results.equil_1.results:z_axis';
-  mapping.expression = 'tcv_eq(''''z_axis'''',''''LIUQE.M'''')';
+  mapping.expression = 'tcv_eq("z_axis","LIUQE.M")';
   %
   % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % extra TCV cases (not necessarily in official data_request name list)
diff --git a/matlab/liuqefortran2liuqematlab.m b/matlab/liuqefortran2liuqematlab.m
index 623d6b520a83785ed10b2fa8b1ce7211ff4883a3..4a337a814222c98ef9707ab403d4547d994bd465 100644
--- a/matlab/liuqefortran2liuqematlab.m
+++ b/matlab/liuqefortran2liuqematlab.m
@@ -8,7 +8,7 @@ function liuqe_fortran_matlab = liuqefortran2liuqematlab(varargin);
 % liuqe_fortran_matlab{:,2} give liuqe matlab names
 %
 % Note that you can have multiple entries (for liuqe fortran since it has less nodes) but the first one should be the closest match (see r_contour/r_rho/r_surf)
-% 
+%
 % varargin: If absent, then return full table
 % varargin{1}: 'node_name_to_convert'
 % varargin{2}: origin of node_name given=varargin{1}: 0: fortran, 1(default): matlab
@@ -27,6 +27,7 @@ function liuqe_fortran_matlab = liuqefortran2liuqematlab(varargin);
 liuqe_fortran_matlab_table = [ ...
     {'l_i'}           , {'l_i_3'}     ; ...
     {'i_p'}           , {'i_pl'}  ; ...
+    {'i_pol_fit'}     , {'i_pol'}  ; ...
     {'pprime_psi'}    , {'pprime_rho'} ; ... % warning on different x-mesh
     {'surface_flux'}  , {'psi_surf'}  ; ...
     {'q_zero'}        , {'q_axis'}    ; ...