add powers and q_rho, changed names with _rho

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

crpptbx_new/TCV/gdat_tcv.m

@@ -236,6 +236,15 @@ if liuqe_version==2 || liuqe_version==3
   substr_liuqe = ['_' num2str(liuqe_version)];
 end
 
+% special treatment for model shot=-1 or preparation shot >=100'000
+begstr = '';
+if shot==-1 || shot>=100000
+  % requires FBTE
+  liuqe_version = -1;
+  begstr = 'tcv_eq( "';
+  substr_liuqe = '", "FBTE" )';
+end
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % Specifications on how to get the data provided in tcv_requests_mapping
@@ -244,8 +253,12 @@ gdat_data.label = mapping_for_tcv.label;
 
 ishot=NaN;
 if do_mdsopen_mdsclose
-  mdsdefaultserver tcv1.epfl.ch; % should be in tcv general path, but set-it in the meantime...
-  ishot = mdsopen(shot); % if ishot equal to shot, then mdsclose at the end
+  % mdsdefaultserver tcv1.epfl.ch; % should be in tcv general path, but set-it in the meantime...
+  if liuqe_version==-1
+    ishot = mdsopen('pcs', shot);
+  else
+    ishot = mdsopen(shot); % if ishot equal to shot, then mdsclose at the end
+  end
   if ishot~=shot
     warning(['cannot open shot= ' num2str(shot)])
     return
@@ -273,7 +286,14 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi')
       return
     end
   else
-    eval_expr = ['tdi(''' mapping_for_tcv.expression substr_tdi ''');'];
+    if liuqe_version==-1
+      if strcmp(lower(mapping_for_tcv.expression(1:8)),'\results')
+        mapping_for_tcv_expression_eff = mapping_for_tcv.expression(11:end);
+      end
+      eval_expr = ['tdi(''' begstr mapping_for_tcv_expression_eff substr_liuqe ''');']
+    else
+      eval_expr = ['tdi(''' mapping_for_tcv.expression substr_tdi ''');'];
+    end
     aatmp=eval(eval_expr);
   end
   if isempty(aatmp.data) || isempty(aatmp.dim) % || ischar(aatmp.data) (to add?)
@@ -403,7 +423,6 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     % B0 at R0=0.88
     R0EXP=0.88;
     if liuqe_version==-1
-      mdsopen( 'pcs', shot); %synthetic shot generated with FBT and MGAMS.
       nodenameeff = 'tcv_eq("BZERO","FBTE")';
       tracetdi=tdi(nodenameeff);
       gdat_data.data = tracetdi.data;
@@ -455,7 +474,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 = beta.dimunits;
+    gdat_data.dimunits = 's';
     
    case {'cxrs'}
     %not yet finished, just started
@@ -697,6 +716,65 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       gdat_data.error_bar = 1.6022e-19 .* (gdat_data.ne.data .* gdat_data.te.error_bar ...
           + 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
+    nodenameeff='\results::toray.input:p_gyro';
+    tracetdi=tdi(nodenameeff);
+    gdat_data.ec.data = tracetdi.data*1e3; % at this stage p_gyro is in kW'
+    gdat_data.ec.units = 'W';
+    gdat_data.ec.dim=tracetdi.dim;
+    gdat_data.ec.dimunits=tracetdi.dimunits;
+    gdat_data.ec.t=tracetdi.dim{1};
+    gdat_data.ec.x=tracetdi.dim{2};
+    gdat_data.ec.data_fullpath=[nodenameeff];
+    gdat_data.ec.label='P_{EC}';
+    % set ec time as reference
+    gdat_data.t = gdat_data.ec.t;
+    gdat_data.dim{1} = gdat_data.t;
+    gdat_data.dimunits{1} = 's';
+    gdat_data.units = 'W';
+    
+    % get ohmic power simply from vloop*Ip (minus sign for TCV)
+    ip=gdat([],'ip');
+    vloop=gdat([],'vloop');
+    tension = -1e5;
+    vloop_smooth=interpos(vloop.t,vloop.data,gdat_data.t,tension);
+    ip_t = interp1(ip.t,ip.data,gdat_data.t);
+    gdat_data.ohm.data = -vloop_smooth.*ip_t;
+    gdat_data.ohm.units = 'W';
+    gdat_data.ohm.dim=gdat_data.dim;
+    gdat_data.ohm.dimunits=gdat_data.dimunits;
+    gdat_data.ohm.t=gdat_data.t;
+    gdat_data.ohm.x=[];
+    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;
+    gdat_data.data_fullpath=['tot power from EC and ohm'];
+    gdat_data.label = 'P_{tot}=P_{EC}+P_{ohm}';
+
+   case {'q_rho'}
+    % q profile on psi from liuqe
+    nodenameeff='\results::q_psi';
+    if liuqe_version==-1
+      nodenameeff=[begstr 'q_psi' substr_liuqe];
+    end
+    tracetdi=tdi(nodenameeff);
+    gdat_data.data = tracetdi.data;
+    gdat_data.dim = tracetdi.dim;
+    gdat_data.t = gdat_data.dim{2};
+    gdat_data.data_fullpath=[nodenameeff ' on rhopol'];
+    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.dimunits{1} = '';
+    gdat_data.dimunits{2} = 's';
+    gdat_data.units = '';
+    gdat_data.request_description = 'q(rhopol\_norm)';
+
    otherwise
     warning(['switchcase= ' data_request_eff ' not known in gdat_tcv'])
     error_status=901;

crpptbx_new/TCV/tcv_requests_mapping.m

@@ -131,6 +131,7 @@ switch lower(data_request)
   mapping.timedim = 2;
   mapping.method = 'switchcase'; % especially since might have option fit, etc
  case 'powers'
+  mapping.timedim = 1;
   mapping.label = 'various powers';
   mapping.method = 'switchcase';
  case 'q0'