From 623f6d84dec46a97f6edc43c35490a62ad5b054a Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Tue, 29 Sep 2015 13:59:39 +0000
Subject: [PATCH] added help for parameters, remove _abs and added _raw for
 thomson ne, added cxrs, fit for conf as default

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@5070 d63d8f72-b253-0410-a779-e742ad2e26cf
---
 crpptbx/TCV/gdat_tcv.m | 261 +++++++++++++++++++++--------------------
 crpptbx/gdat.m         |   2 +
 2 files changed, 136 insertions(+), 127 deletions(-)

diff --git a/crpptbx/TCV/gdat_tcv.m b/crpptbx/TCV/gdat_tcv.m
index ee7f6b09..c0381d59 100644
--- a/crpptbx/TCV/gdat_tcv.m
+++ b/crpptbx/TCV/gdat_tcv.m
@@ -33,9 +33,11 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(shot,data_req
 % gdat_data.machine: machine providing the data
 % gdat_data.gdat_request: keyword for gdat if relevant
 % gdat_data.data_fullpath: full path to the data node if known and relevant, or relevant expression called if relevant
-% gdat_data.gdat_params: copy gdat_params for completeness
+% gdat_data.gdat_params: copy gdat_params for completeness (gdat_params contains a .help structure detailing each parameter)
 % gdat_data.xxx: any other relevant information
 %
+% gdat_params contains the options relevant for the called data_request. It also contains a help structure for each option
+% eg.: param1 in gdat_params.param1 and help in gdat_params.help.param1
 %
 % Examples:
 % (should add working examples for various machines (provides working shot numbers for each machine...))
@@ -105,6 +107,9 @@ data_request_eff = '';
 if nargin>=2 && ischar(data_request); data_request = lower(data_request); end
 
 gdat_data.gdat_request = data_request_names_all; % so if return early gives list of possible request names
+gdat_data.gdat_params.help = tcv_help_parameters(fieldnames(gdat_data.gdat_params));
+gdat_params = gdat_data.gdat_params;
+
 % no inputs
 if nargin==0
   % return defaults and list of keywords
@@ -427,6 +432,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % First the request names valid for "all" machines:
     %
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'a_minor','rgeom'}
     % compute average minor or major radius (on z=zaxis normally)
     nodenameeff=['\results::r_max_psi' substr_liuqe];
@@ -454,6 +460,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
     gdat_data.dimunits = rmaxpsi.dimunits(2);
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'zgeom'}
     % compute average minor or major radius (on z=zaxis normally)
     nodenameeff=['\results::z_contour' substr_liuqe];
@@ -472,6 +479,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       gdat_data.units = zcontour.units;
     end
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'b0'}
     % B0 at R0=0.88
     R0EXP=0.88;
@@ -497,6 +505,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.dimunits = tracetdi.dimunits;
     gdat_data.request_description = ['vacuum magnetic field at R0=' num2str(R0EXP) 'm; COCOS=17'];
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'betan'}
     % 100*beta / |Ip[MA] * B0[T]| * a[m]
     % get B0 from gdat_tcv, without re-opening the shot and using the same parameters except data_request
@@ -529,6 +538,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.units = '';
     gdat_data.dimunits{1} = 's';
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'cxrs'}
     %not yet finished, just started
     % load typical data from cxrs, Ti, ni, vtori and vpoli (if available), as well as zeff from cxrs
@@ -618,6 +628,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.dimunits{1} = '';
     gdat_data.dimunits{2} = 's';
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'eqdsk'}
     %
     time=1.; % default time
@@ -690,6 +701,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.request_description = ['data=psi, x=(R,Z), eqdsk contains eqdsk structure with which ' ...
                     'plot_eqdsk, write_eqdsk, read_eqdsk can be used'];
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'mhd'}
     % load n=1, 2 and 3 Bdot from magnetic measurements
     n1=tdi('abs(mhdmode("LFS",1,1))');
@@ -709,115 +721,36 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       gdat_data.request_description = 'delta_Bdot from magnetic probes to get n=1, 2 and 3';
     end
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ne','te'}
     % ne or Te from Thomson data on raw z mesh vs (z,t)
-    edge_str_ = '';
-    edge_str_dot = '';
-    if isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
-          gdat_data.gdat_params.edge>0
-      edge_str_ = '_edge';
-      edge_str_dot = '.edge';
-    else
+    if ~(isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
+         gdat_data.gdat_params.edge>0)
       gdat_data.gdat_params.edge = 0;
     end
-    nodenameeff=['\results::thomson' edge_str_dot ':' data_request_eff];
-    tracetdi=tdi(nodenameeff);
-    tracestd=tdi(['\results::thomson' edge_str_dot ':' data_request_eff ':error_bar']);
-    gdat_data.data=tracetdi.data'; % Thomson data as (t,z)
-    gdat_data.error_bar=tracestd.data';
-    gdat_data.data_fullpath=[nodenameeff];
-    % add correct dimensions
-    try
-      time=mdsdata('\results::thomson:times');
-    catch
-      if (nverbose>=1) && shot<100000
-        warning('Problems with \results::thomson:times')
-        disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
-      end
-      return
-    end
-    if isempty(time) || ischar(time)
-      thomsontimes=time;
-      if (nverbose>=1) && shot<100000
-        warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
-        disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
-      end
-      return
-    end
-    if strcmp(data_request_eff(1:2),'ne')
-      tracefirrat_data = get_fir_thom_rat_data(shot,['thomson' edge_str_],time);
-      gdat_data.data_abs = gdat_data.data * diag(tracefirrat_data);
-      gdat_data.error_bar_abs = gdat_data.error_bar * diag(tracefirrat_data);
-      gdat_data.firrat=tracefirrat_data;
-      gdat_data.data_fullpath=[gdat_data.data_fullpath ' ; _abs includes *firrat'];
-    end
-    z=mdsdata('\diagz::thomson_set_up:vertical_pos');
-    gdat_data.dim=[{z};{time}];
-    gdat_data.dimunits=[{'Z [m]'} ; {'time [s]'}];
-    gdat_data.x=z;
-    gdat_data.t=time;
-    % isfield does not work since tracetdi is not a 'struct' but a tdi object, thus use fieldnames
-    if any(strcmp(fieldnames(tracetdi),'units'))
-      gdat_data.units=tracetdi.units;
-    end
+    [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.edge,nverbose);
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ne_rho', 'te_rho', 'nete_rho'}
-    try
-      time=mdsdata('\results::thomson:times');
-    catch
-      if (nverbose>=1) && shot<100000
-        warning('Problems with \results::thomson:times')
-        warning(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
-      end
-      return
-    end
-    if isempty(time) || ischar(time)
-      thomsontimes=time;
-      if (nverbose>=1) && shot<100000
-        warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
-        disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
-      end
-      return
-    end
+    % if nete_rho, do first ne, then Te later (so fir stuff already done)
     zshift = 0.;
     if isfield(gdat_data.gdat_params,'zshift') && ~isempty(gdat_data.gdat_params.zshift)
       zshift = gdat_data.gdat_params.zshift;
     else
       gdat_data.gdat_params.zshift = zshift;
     end
+    if ~(isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
+         gdat_data.gdat_params.edge>0)
+      gdat_data.gdat_params.edge = 0;
+    end
+    [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.edge,nverbose);
+    % construct rho mesh
     edge_str_ = '';
     edge_str_dot = '';
-    if isfield(gdat_data.gdat_params,'edge') && ~isempty(gdat_data.gdat_params.edge) && ...
-          gdat_data.gdat_params.edge>0
+    if gdat_data.gdat_params.edge
       edge_str_ = '_edge';
       edge_str_dot = '.edge';
-    else
-      gdat_data.gdat_params.edge = 0;
-    end
-    % if nete_rho, do first ne, then Te later (so fir stuff already done)
-    if strcmp(data_request_eff,'ne_rho')  || strcmp(data_request_eff,'nete_rho')
-      nodenameeff=['\results::thomson' edge_str_dot ':ne'];
-      tracetdi=tdi(nodenameeff);
-      nodenameeff=['\results::thomson' edge_str_dot ':ne; error_bar ; fir_thom_rat; (ne,std)*fir_thom_rat'];
-      tracestd=tdi(['\results::thomson'  edge_str_dot ':ne:error_bar']);
-      tracefirrat_data = get_fir_thom_rat_data(shot,['thomson' edge_str_],time);
-    else
-      nodenameeff=['\results::thomson' edge_str_dot ':te'];
-      tracetdi=tdi(nodenameeff);
-      nodenameeff=['\results::thomson' edge_str_dot ':te; error_bar'];
-      tracestd=tdi(['\results::thomson' edge_str_dot ':te:error_bar']);
     end
-    gdat_data.data=tracetdi.data'; % Thomson data as (t,z)
-    gdat_data.error_bar=tracestd.data';
-    gdat_data.data_fullpath=nodenameeff;
-    if strcmp(data_request_eff,'ne_rho')  || strcmp(data_request_eff,'nete_rho')
-      gdat_data.firrat=tracefirrat_data;
-      gdat_data.data_abs=gdat_data.data*diag(tracefirrat_data);
-      gdat_data.error_bar_abs=gdat_data.error_bar*diag(tracefirrat_data);
-      gdat_data.data_fullpath=[gdat_data.data_fullpath ' ; _abs includes *firrat'];
-    end
-    % add correct dimensions
-    % construct rho mesh
     psi_max=tdi(['\results::thomson' edge_str_dot ':psi_max' substr_liuqe]);
     psiscatvol=tdi(['\results::thomson' edge_str_dot ':psiscatvol' substr_liuqe]);
     if abs(zshift)>1e-5
@@ -839,6 +772,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
           disp(' bad time dimension for zshift')
           disp(['it should be 1 or ' num2str(length(psiscatvol.dim{1})) ' or ' num2str(length(psitdi.dim{3}))])
       end
+max(zeffshift)
       for it=1:length(psiscatvol.dim{1})
         itpsitdi=iround(psitdi.dim{3},psiscatvol.dim{1}(it));
         psirz=psitdi.data(:,:,itpsitdi);
@@ -853,36 +787,42 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     else
       rho=NaN;
     end
-    gdat_data.dim=[{rho};{time}];
-    gdat_data.dimunits=[{'sqrt(psi)'} ; {'time [s]'}];
+    gdat_data.dim{1}=rho;
+    gdat_data.dimunits=[{'sqrt(psi_norm)'} ; {'time [s]'}];
     gdat_data.x=rho;
-    gdat_data.t=time;
-    if any(strcmp(fieldnames(tracetdi),'units'))
-      gdat_data.units=tracetdi.units;
-    end
     %%%%%%%%%%% add fitted profiles if 'fit'>=1
+    default_fit_type = 'conf';
     if isfield(gdat_data.gdat_params,'fit') && ~isempty(gdat_data.gdat_params.fit) && ...
           gdat_data.gdat_params.fit>0
       % default is from proffit:avg_time
-      def_proffit = '\results::proffit.avg_time';
-      if isfield(gdat_data.gdat_params,'fit_type') && ~isempty(gdat_data.gdat_params.fit_type)
-        if strcmp(gdat_data.gdat_params.fit_type,'local')
-          def_proffit = '\results::proffit.local_time';
-        else
-          gdat_data.gdat_params.fit_type = 'avg';
-        end
-      else
-        gdat_data.gdat_params.fit_type = 'avg';
+      if ~(isfield(gdat_data.gdat_params,'fit_type') && ~isempty(gdat_data.gdat_params.fit_type)) || ~any(strcmp(gdat_data.gdat_params.fit_type,{'local', 'avg', 'conf'}))
+        gdat_data.gdat_params.fit_type = default_fit_type;
       end
-      if strcmp(data_request_eff(1:2),'ne')
-        nodenameeff = [def_proffit ':neft_abs']; % do first ne if nete asked for
-      elseif strcmp(data_request_eff(1:2),'te')
-        nodenameeff = [def_proffit ':teft'];
+      switch gdat_data.gdat_params.fit_type
+       case 'avg'
+        def_proffit = '\results::proffit.avg_time:';
+       case 'local'
+        def_proffit = '\results::proffit.local_time:';
+       case 'conf'
+        def_proffit = '\results::conf:';
+       otherwise
+        disp('should not be in switch gdat_data.gdat_params.fit_type')
+        disp('ask olivier.sauter@epfl.ch')
+        return
+      end
+      if strcmp(gdat_data.gdat_params.fit_type,'conf')
+        nodenameeff = [def_proffit data_request_eff(1:2)];
       else
-        disp(['should not be here: data_request_eff, data_request_eff(1:2)= ',data_request_eff, data_request_eff(1:2)]);
+        if strcmp(data_request_eff(1:2),'ne')
+          nodenameeff = [def_proffit 'neft_abs']; % do first ne if nete asked for
+        elseif strcmp(data_request_eff(1:2),'te')
+          nodenameeff = [def_proffit 'teft'];
+        else
+          disp(['should not be here: data_request_eff, data_request_eff(1:2)= ',data_request_eff, data_request_eff(1:2)]);
+        end
       end
       if isfield(gdat_data.gdat_params,'trialindx') && ~isempty(gdat_data.gdat_params.trialindx) && ...
-          gdat_data.gdat_params.trialindx>=0
+            gdat_data.gdat_params.trialindx>=0
         nodenameeff=[nodenameeff ':trial'];
         trialindx = gdat_data.gdat_params.trialindx;
       else
@@ -905,8 +845,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       gdat_data.fit.t=tracetdi.dim{2};
       if mapping_for_tcv.timedim~=2 | mapping_for_tcv.gdat_timedim~=2
         disp(['unexpected timedim in fit: data_request_eff= ' data_request_eff ...
-             ', mapping_for_tcv.timedim= ' mapping_for_tcv.timedim ...
-             ', mapping_for_tcv.gdat_timedim= ' mapping_for_tcv.gdat_timedim]);
+              ', mapping_for_tcv.timedim= ' mapping_for_tcv.timedim ...
+              ', mapping_for_tcv.gdat_timedim= ' mapping_for_tcv.gdat_timedim]);
       end
       gdat_data.dim=tracetdi.dim(1:2);
       gdat_data.dimunits=tracetdi.dimunits(1:2);
@@ -918,17 +858,21 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       if strcmp(data_request_eff(1:4),'nete')
         gdat_data.fit.ne.data = gdat_data.fit.data;
         gdat_data.fit.ne.units = gdat_data.fit.units;
-        nodenameeff = [def_proffit ':teft'];
+        if strcmp(gdat_data.gdat_params.fit_type,'conf')
+          nodenameeff = [def_proffit 'te'];
+        else
+          nodenameeff = [def_proffit 'teft'];
+        end
         if ~isempty(trialindx); nodenameeff=[nodenameeff ':trial']; end
         tracetdi=tdi(nodenameeff);
         if isempty(trialindx)
           gdat_data.fit.te.data = tracetdi.data;
         else
-        if ~isempty(tracetdi.data) && size(tracetdi.data,3)>=trialindx+1
-          gdat_data.fit.te.data = tracetdi.data(:,:,trialindx+1);
-        else
-          return
-        end
+          if ~isempty(tracetdi.data) && size(tracetdi.data,3)>=trialindx+1
+            gdat_data.fit.te.data = tracetdi.data(:,:,trialindx+1);
+          else
+            return
+          end
         end
         if any(strcmp(fieldnames(tracetdi),'units'))
           gdat_data.fit.te.units=tracetdi.units;
@@ -940,15 +884,16 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       end
     else
       gdat_data.gdat_params.fit = 0;
+      gdat_data.gdat_params.fit_type = default_fit_type;
     end
     %%%%%%%%%%%
     % if nete_rho, copy data as .ne, get .te and put pe=e ne Te in data:
     if strcmp(data_request_eff(1:4),'nete')
-      gdat_data.ne.data = gdat_data.data_abs;
-      gdat_data.ne.error_bar = gdat_data.error_bar_abs;
+      gdat_data.ne.data = gdat_data.data;
+      gdat_data.ne.error_bar = gdat_data.error_bar;
       gdat_data.ne.firrat=gdat_data.firrat;
       gdat_data.ne.units = 'm^{-3}';
-      gdat_data = rmfield(gdat_data,{'firrat','data_abs','error_bar_abs'});
+      gdat_data = rmfield(gdat_data,{'firrat','data_raw','error_bar_raw'});
       %
       nodenameeff=['\results::thomson' edge_str_dot ':te'];
       tracetdi=tdi(nodenameeff);
@@ -965,6 +910,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
           + gdat_data.te.data .* gdat_data.ne.error_bar);
     end
 
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'powers'}
     % note: same time array for all main, ec, ohm, nbi, ...
     % At this stage fill just ech, later add nbi
@@ -1008,6 +954,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.data_fullpath=['tot power from EC and ohm'];
     gdat_data.label = 'P_{ohm};P_{EC};P_{tot}';
 
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'q_rho'}
     % q profile on psi from liuqe
     nodenameeff=['\results::q_psi' substr_liuqe];
@@ -1028,6 +975,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.units = '';
     gdat_data.request_description = 'q(rhopol\_norm)';
 
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'psi_edge'}
     % psi at edge, 0 by construction in Liuqe, thus not given
     nodenameeff=['\results::psi_axis' substr_liuqe];
@@ -1043,6 +991,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.units = tracetdi.units;
     gdat_data.request_description = '0 since LIUQE construct psi to be zero at LCFS';
 
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'rhotor_edge','rhotor'}
     % Phi(LCFS) = int(Bphi dSphi), can use Eq.(11) of "Tokamak coordinate conventions: COCOS" paper:
     % O. Sauter, S.Yu. Medvedev, Comput. Phys. Commun. 184 (2013) 293–302
@@ -1096,6 +1045,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       gdat_data.b0 = b0tpsi(it);
     end
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'rhovol','volume_rho','volume'}
     % volume_rho = vol(rho); volume = vol(LCFS) = vol(rho=1);
     % rhovol = sqrt(vol(rho)/vol(rho=1));
@@ -1133,9 +1083,11 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
       gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'sxr'}
     % sxr from Xtomo by default or dmpx if 'camera','dmpx' is provided
 
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'profnerho','profterho'}
     % for backward compatibility but corresponds to ne_rho with param.fit_type='auto' (TCV special)
     % 
@@ -1201,7 +1153,9 @@ end
 
 if ishot==shot; mdsclose; end
 
+gdat_data.gdat_params.help = tcv_help_parameters(fieldnames(gdat_data.gdat_params));
 gdat_data.mapping_for.tcv = mapping_for_tcv;
+gdat_params = gdat_data.gdat_params;
 error_status=0;
 
 return
@@ -1262,8 +1216,61 @@ else
   disp('bad input in get_fir_thom_rat_data')
   return
 end  
-  
-if ~isempty(tracefirrat.data) || ischar(tracefirrat.data)
+
+if ~isempty(tracefirrat.data) && ~ischar(tracefirrat.data) && any(~isnan(tracefirrat.data)) ...
+      && ~isempty(tracefirrat.dim) && ~isempty(tracefirrat.dim{1})
   firthomratio = interp1(tracefirrat.dim{1},tracefirrat.data,timebase);
 end
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,doedge,nverbose);
+%
+try
+  time=mdsdata('\results::thomson:times');
+catch
+  if (nverbose>=1) && shot<100000
+    warning('Problems with \results::thomson:times')
+    warning(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+  end
+  return
+end
+if isempty(time) || ischar(time)
+  thomsontimes=time;
+  if (nverbose>=1) && shot<100000
+    warning('!!!!!!!!!!!!!!!!!!!!!!!!!\results::thomson:times  is empty? Check')
+    disp(['!!!!!!!!!!!!!!!!!!!!!!!!! cannot continue with ' data_request_eff])
+  end
+  return
+end
+edge_str_ = '';
+edge_str_dot = '';
+if doedge
+  edge_str_ = '_edge';
+  edge_str_dot = '.edge';
+end
+
+nodenameeff = ['\results::thomson' edge_str_dot ':' data_request_eff(1:2)];
+tracetdi=tdi(nodenameeff);
+gdat_data.data=tracetdi.data'; % Thomson data as (t,z)
+tracestd=tdi(['\results::thomson'  edge_str_dot ':' data_request_eff(1:2) ':error_bar']);
+gdat_data.error_bar=tracestd.data';
+gdat_data.data_fullpath=[nodenameeff '; error_bar'];
+% add fir if ne requested
+if strcmp(data_request_eff(1:2),'ne')
+  tracefirrat_data = get_fir_thom_rat_data(shot,['thomson' edge_str_],time);
+  gdat_data.firrat=tracefirrat_data;
+  gdat_data.data_raw = gdat_data.data;
+  gdat_data.data = gdat_data.data_raw * diag(tracefirrat_data);
+  gdat_data.error_bar_raw = gdat_data.error_bar;
+  gdat_data.error_bar = gdat_data.error_bar_raw * diag(tracefirrat_data);
+  gdat_data.data_fullpath=[gdat_data.data_fullpath '; fir_thom_rat ; _raw without firrat'];
+end
+z=mdsdata('\diagz::thomson_set_up:vertical_pos');
+gdat_data.dim=[{z};{time}];
+gdat_data.dimunits=[{'Z [m]'} ; {'time [s]'}];
+gdat_data.x=z;
+gdat_data.t=time;
+% isfield does not work since tracetdi is not a 'struct' but a tdi object, thus use fieldnames
+if any(strcmp(fieldnames(tracetdi),'units'))
+  gdat_data.units=tracetdi.units;
+end
diff --git a/crpptbx/gdat.m b/crpptbx/gdat.m
index 5b4804f9..5d40f667 100644
--- a/crpptbx/gdat.m
+++ b/crpptbx/gdat.m
@@ -37,6 +37,8 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request
 % gdat_data.gdat_params: copy gdat_params for completeness
 % gdat_data.xxx: any other relevant information
 %
+% gdat_params contains the options relevant for the called data_request. It also contains a help structure for each option
+% eg.: param1 in gdat_params.param1 and help in gdat_params.help.param1
 %
 % Examples:
 % (should add working examples for various machines (provides working shot numbers for each machine...))
-- 
GitLab