add several cases

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

crpptbx_new/TCV/gdat_tcv.m

@@ -287,6 +287,7 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi')
     end
   else
     if liuqe_version==-1
+      mapping_for_tcv_expression_eff = mapping_for_tcv.expression;
       if strcmp(lower(mapping_for_tcv.expression(1:8)),'\results')
         mapping_for_tcv_expression_eff = mapping_for_tcv.expression(11:end);
       end
@@ -401,7 +402,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % First the request names valid for "all" machines:
     %
-   case {'a_minor'}
+   case {'a_minor','rgeom'}
+    % compute average minor or major radius (on z=zaxis normally)
     nodenameeff=['\results::r_max_psi' substr_liuqe];
     rmaxpsi=tdi(nodenameeff);
     nodenameeff2=['\results::r_min_psi' substr_liuqe];
@@ -410,8 +412,16 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     if ~isempty(ij); rmaxpsi.data(ij)=NaN; end
     ij=find(rminpsi.data<0.5 | rminpsi.data>1.2);
     if ~isempty(ij); rminpsi.data(ij)=NaN; end
-    gdat_data.data=0.5.*(rmaxpsi.data(end,:) - rminpsi.data(end,:));
-    gdat_data.data_fullpath=[nodenameeff ' - ' nodenameeff2 ' /2'];
+    if strcmp(data_request_eff,'a_minor')
+      gdat_data.data=0.5.*(rmaxpsi.data(end,:) - rminpsi.data(end,:));
+      gdat_data.data_fullpath=[nodenameeff ' - ' nodenameeff2 ' /2'];
+    elseif strcmp(data_request_eff,'rgeom')
+      gdat_data.data=0.5.*(rmaxpsi.data(end,:) + rminpsi.data(end,:));
+      gdat_data.data_fullpath=[nodenameeff ' + ' nodenameeff2 ' /2'];
+    else
+      disp(['should not be in this case with data_request_eff = ' data_request_eff])
+      return
+    end
     gdat_data.dim = rmaxpsi.dim(2);    
     gdat_data.t = gdat_data.dim{1};
     if any(strcmp(fieldnames(rmaxpsi),'units'))
@@ -419,6 +429,24 @@ 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];
+    zcontour=tdi(nodenameeff);
+    if strcmp(data_request_eff,'zgeom')
+      gdat_data.data=0.5.*(max(zcontour.data,[],1) + min(zcontour.data,[],1));
+      gdat_data.data_fullpath=['(max+min)/2 of ' nodenameeff];
+      gdat_data.dim{1} = zcontour.dim{2};
+      gdat_data.dimunits{1} = zcontour.dimunits{2};
+    else
+      disp(['should not be in this case with data_request_eff = ' data_request_eff])
+      return
+    end
+    gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
+    if any(strcmp(fieldnames(zcontour),'units'))
+      gdat_data.units = zcontour.units;
+    end
+    
    case {'b0'}
     % B0 at R0=0.88
     R0EXP=0.88;
@@ -474,7 +502,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.data = 100.*beta.data ./ abs(ip_t).*1.e6 .* abs(b0_t) .* a_minor_t;
     gdat_data.data_fullpath='100*beta/ip*1e6*b0*a_minor, each from gdat_tcv';
     gdat_data.units = '';
-    gdat_data.dimunits = 's';
+    gdat_data.dimunits{1} = 's';
     
    case {'cxrs'}
     %not yet finished, just started
@@ -751,14 +779,18 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.ohm.data_fullpath=['-vloop(tens=' num2str(tension,'%.0e') ')*ip, from gdat'];
     gdat_data.ohm.label='P_{OHM}';
     
-    % total power from all above
-    gdat_data.data = gdat_data.ec.data(:,10) + gdat_data.ohm.data;
+    % total power from each and total
+    gdat_data.data(:,1) = gdat_data.ohm.data;
+    gdat_data.data(:,2) = gdat_data.ec.data(:,10);
+    gdat_data.data(:,3) = gdat_data.ec.data(:,10) + gdat_data.ohm.data;
+    gdat_data.dim{2} = [1:3];
+    gdat_data.dimunits{2} = 'Pohm;Pec;Ptot';
     gdat_data.data_fullpath=['tot power from EC and ohm'];
-    gdat_data.label = 'P_{tot}=P_{EC}+P_{ohm}';
+    gdat_data.label = 'P_{ohm};P_{EC};P_{tot}';
 
    case {'q_rho'}
     % q profile on psi from liuqe
-    nodenameeff='\results::q_psi';
+    nodenameeff=['\results::q_psi' substr_liuqe];
     if liuqe_version==-1
       nodenameeff=[begstr 'q_psi' substr_liuqe];
     end
@@ -770,11 +802,117 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     rhopol_eff = ones(size(tracetdi.dim{1}));
     rhopol_eff(:) = sqrt(linspace(0,1,length(tracetdi.dim{1})));
     gdat_data.dim{1} = rhopol_eff;
+    gdat_data.x = gdat_data.dim{1};
     gdat_data.dimunits{1} = '';
     gdat_data.dimunits{2} = 's';
     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];
+    if liuqe_version==-1
+      nodenameeff=[begstr 'q_psi' substr_liuqe];
+    end
+    tracetdi=tdi(nodenameeff);
+    gdat_data.data = tracetdi.data.*0;
+    gdat_data.dim = tracetdi.dim;
+    gdat_data.t = gdat_data.dim{1};
+    gdat_data.data_fullpath=[' zero '];
+    gdat_data.dimunits = tracetdi.dimunits;
+    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
+    % since cocos=17 for LIUQE we get:
+    % q = -dPhi/dpsi => Phi = - int(q*dpsi) which should always have the sign of B0
+    params_eff = gdat_data.gdat_params;
+    params_eff.data_request='q_rho'; 
+    q_rho=gdat_tcv([],params_eff);
+    params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE
+    psi_axis=gdat_tcv([],params_eff);
+    params_eff.data_request='b0'; % psi_edge=0 with LIUQE
+    b0=gdat_tcv([],params_eff);
+    b0tpsi = interp1(b0.t,b0.data,psi_axis.t); %q_rho on same time base as psi_axis
+    if isempty(psi_axis.data) || isempty(psi_axis.dim) || isempty(q_rho.data) || isempty(q_rho.dim)
+      warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff])
+      return
+    end
+    rhoequal = linspace(0,1,length(q_rho.dim{1}));
+    if strcmp(data_request,'rhotor_edge')
+      gdat_data.data = psi_axis.data; % to have the dimensions correct
+      gdat_data.dim = psi_axis.dim;
+      gdat_data.t = gdat_data.dim{1};
+      gdat_data.data_fullpath='phi from q_rho, psi_axis and integral(-q dpsi)';
+      gdat_data.units = 'T m^2';
+      gdat_data.dimunits{1} = 's';
+    elseif strcmp(data_request,'rhotor')
+      gdat_data.data = q_rho.data; % to have the dimensions correct
+      gdat_data.dim{1} = ones(size(q_rho.dim{1}));
+      gdat_data.dim{1}(:) = rhoequal;
+      gdat_data.dim{2} = q_rho.dim{2};
+      gdat_data.t = gdat_data.dim{2};
+      gdat_data.data_fullpath='sqrt(phitor/phitor_edge), rhotor_edge=sqrt(phitor/B0/pi)';
+      gdat_data.units = '';
+      gdat_data.dimunits{1} = 'rhopol\_norm';
+      gdat_data.dimunits{2} = 's';
+    end
+    for it=1:length(psi_axis.data)
+      ij=find(~isnan(q_rho.data(:,it)));
+      if ~isempty(ij)
+        [qfit,~,~,phi]=interpos(q_rho.x(ij).^2,q_rho.data(ij,it),rhoequal.^2);
+        dataeff = sqrt(phi .* psi_axis.data(it) ./ b0tpsi(it) ./ pi) ; % Delta_psi = -psi_axis
+      else
+        dataeff = NaN;
+      end
+      if strcmp(data_request,'rhotor_edge')
+        gdat_data.data(it) = dataeff(end);
+      elseif strcmp(data_request,'rhotor')
+        gdat_data.data(:,it) = dataeff./dataeff(end);
+        gdat_data.rhotor_edge(it) = dataeff(end);
+      end
+      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));
+    nodenameeff='\results::psitbx:vol';
+    if liuqe_version==-1
+      nodenameeff=[begstr 'vol' substr_liuqe];
+    end
+    tracetdi=tdi(nodenameeff);
+    if isempty(tracetdi.data) || isempty(tracetdi.dim) 
+      return
+    end
+    gdat_data.units = tracetdi.units;
+    if strcmp(data_request,'volume')
+      gdat_data.data = tracetdi.data(end,:);
+      gdat_data.dim{1} = tracetdi.dim{2};
+      gdat_data.data_fullpath=['\results::psitbx:vol(end,:)'];
+      gdat_data.dimunits{1} = tracetdi.dimunits{2};
+      gdat_data.request_description = 'volume(LCFS)=volume(rhopol=1)';
+    else
+      gdat_data.data = tracetdi.data;
+      gdat_data.dim = tracetdi.dim;
+      gdat_data.dimunits = tracetdi.dimunits;
+      if strcmp(data_request,'volume_rho')
+        gdat_data.data_fullpath=['\results::psitbx:vol'];
+        gdat_data.request_description = 'volume(rho)';
+      elseif strcmp(data_request,'rhovol')
+        gdat_data.volume_edge = gdat_data.data(end,:);
+        gdat_data.data = sqrt(gdat_data.data./repmat(reshape(gdat_data.volume_edge,1,size(gdat_data.data,2)),size(gdat_data.data,1),1));
+        gdat_data.data_fullpath='sqrt(\results::psitbx:vol/vol_edge)';
+        gdat_data.request_description = 'sqrt(volume(rho)/volume(edge))';
+      else
+        disp(['should not be here in vol cases with data_request = ' data_request]);
+        return
+      end
+    end
+      gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim};
+    
    otherwise
     warning(['switchcase= ' data_request_eff ' not known in gdat_tcv'])
     error_status=901;

crpptbx_new/TCV/tcv_requests_mapping.m

@@ -134,6 +134,14 @@ switch lower(data_request)
   mapping.timedim = 1;
   mapping.label = 'various powers';
   mapping.method = 'switchcase';
+ case 'psi_axis'
+  mapping.timedim = 1;
+  mapping.method = 'tdiliuqe';
+  mapping.expression = '\results::psi_axis';
+ case 'psi_edge'
+  mapping.timedim = 1;
+  mapping.method = 'switchcase'; % is set to zero, so not in tree nodes
+  mapping.label = 'psi\_edge';
  case 'q0'
   mapping.timedim = 1;
   mapping.method = 'tdiliuqe';
@@ -191,7 +199,7 @@ switch lower(data_request)
   mapping.method = 'switchcase';
  case 'vloop'
   mapping.timedim = 1;
-  mapping.label = '';
+  mapping.label = 'Vloop';
   mapping.method = 'tdi';
   mapping.expression = [{'\magnetics::vloop[*,$1]'} {'001'}];
  case 'volume'

crpptbx_new/test_all_requestnames.m

+
+machine='TCV';
+shot=48836;
+aa=gdat('machine',machine);
+all_request_names = aa.gdat_request
+%break
+istart=1;
+for irequest=istart:length(all_request_names)
+  request=all_request_names{irequest}
+  ab{irequest} = gdat(shot,request,'machine',machine,'doplot',1);
+  pause
+end