From 0e9877a6bb9bc89ad55f39f8591f22f2449f098e Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Fri, 28 Jun 2019 09:03:08 +0000
Subject: [PATCH] add Plh for TCV and JET and other jet stuff
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@12238 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/JET/gdat_jet.m | 55 +++++++++++++++---------------
crpptbx/JET/jet_requests_mapping.m | 25 +++++++++++++-
crpptbx/TCV/tcv_requests_mapping.m | 31 +++++++++++++++++
3 files changed, 83 insertions(+), 28 deletions(-)
diff --git a/crpptbx/JET/gdat_jet.m b/crpptbx/JET/gdat_jet.m
index 5a17ffba..ba994cfe 100644
--- a/crpptbx/JET/gdat_jet.m
+++ b/crpptbx/JET/gdat_jet.m
@@ -2,7 +2,7 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_jet(shot,data_req
%
% function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
%
-% Aim: get data from a given machine using full path or keywords.
+% Aim: get data from a given machine using full path or keywords.
% data_request are and should be case independent (transformed in lower case in the function and outputs)
%
% If no inputs are provided, return the list of available pre-defined data_request in gdat_data and default parameters gdat_params
@@ -41,7 +41,7 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_jet(shot,data_req
%
% Examples:
% (should add working examples for various machines (provides working shot numbers for each machine...))
-%
+%
% [a1,a2]=gdat;
% a2.data_request = 'Ip';
% a3=gdat(51994,a2); % gives input parameters as a structure, allows to call the same for many shots
@@ -209,7 +209,7 @@ if ischar(data_request_eff) || length(data_request_eff)==1
else
ij=[];
end
-if ~isempty(ij);
+if ~isempty(ij);
gdat_data.gdat_request = data_request_names_all{ij};
if isfield(data_request_names.all.(data_request_names_all{ij}),'description') && ~isempty(data_request_names.all.(data_request_names_all{ij}).description)
% copy description of keyword
@@ -310,7 +310,7 @@ if strcmp(mapping_for_jet.method,'signal')
% keyboard
if mapping_for_jet.gdat_timedim>0 && mapping_for_jet.gdat_timedim ~= mapping_for_jet.timedim
% shift timedim to gdat_timedim data(i,j,...itime,k,...) -> data(i,inewtime,j,...,k,...)
- % note that this means that gdat_data.x and gdat_data.t are same and correct,
+ % note that this means that gdat_data.x and gdat_data.t are same and correct,
% only .data, .dim and .dimunits need to be changed
iprev=[1:nbdims];
ij=find(dim_nontim>mapping_for_jet.gdat_timedim-1);
@@ -359,7 +359,7 @@ elseif strcmp(mapping_for_jet.method,'expression')
ijdim=find(strcmp(tmp_fieldnames,'dim')==1);
if ~isempty(ijdim)
nbdims = length(gdat_data.dim);
- if mapping_for_jet.timedim==-1;
+ if mapping_for_jet.timedim==-1;
mapping_for_jet.timedim = nbdims;
if (size(gdat_data.data,nbdims)==1 && nbdims>1); mapping_for_jet.timedim = nbdims-1; end
end
@@ -390,16 +390,16 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
% First the request names valid for "all" machines:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
% if 'fit' option is added: 'fit',1, then the fitted profiles are returned
- %
+ %
% sub_nodes names from CXRS_get_profiles function, lower case when used in gdat
sub_nodes = {'Ti','vTor','vPol','nC','Zeff'}; % first node is also copied into data, choose "default' one
- sub_nodes_out = lower({'Ti','vTor','vPol','ni','Zeff'}); %
+ sub_nodes_out = lower({'Ti','vTor','vPol','ni','Zeff'}); %
sub_nodes_units = {'eV','m/s','m/s','m^{-3}',''}; % first node is also copied into data, choose "default' one
% use A. Karpushov routine to get profiles and then copy the data or the fitted profiles
aa=CXRS_get_profiles; cxrs_params = aa.param;
@@ -448,7 +448,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
cxrs_params.prof.Ti.taus = fit_tension.ti;
cxrs_params.prof.vi.taus = fit_tension.vi;
cxrs_params.prof.nc.taus = fit_tension.ni;
- cxrs_params.prof.zeff.taus = fit_tension.zeff;
+ cxrs_params.prof.zeff.taus = fit_tension.zeff;
cxrs_params.k_plot = cxrs_plot;
cxrs_profiles = CXRS_get_profiles(shot,source,time_interval,cxrs_params);
inb_times = length(cxrs_profiles.Times);
@@ -506,7 +506,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
gdat_data.(sub_eff_out).raw.error_bar_rho(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.dx;
gdat_data.(sub_eff_out).raw.cxrs_system(1:nb_raw_points,it) = cxrs_profiles.(sub_eff){it}.use.sys;
gdat_data.time_interval{it} = cxrs_profiles.(sub_eff){it}.t_lim;
- end
+ end
end
gdat_data.(sub_eff_out).units = sub_nodes_units{i};
if i==1
@@ -525,7 +525,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
end
gdat_data.dimunits{1} = '';
gdat_data.dimunits{2} = 's';
-
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'eqdsk'}
%
@@ -558,7 +558,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
gdat_data.dim{1} = gdat_data.eqdsk{itime}.rmesh;
gdat_data.dim{2} = gdat_data.eqdsk{itime}.zmesh;
else
- % for several times, use array of structure for eqdsks,
+ % for several times, use array of structure for eqdsks,
% cannot use it for psi(R,Z) in .data and .x since R, Z might be different at different times,
% so project psi(R,Z) on Rmesh, Zmesh of 1st time
xx=repmat(reshape(gdat_data.dim{1},length(gdat_data.dim{1}),1),1,size(gdat_data.eqdsk{itime}.psi,2));
@@ -583,7 +583,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
'plot_eqdsk, write_eqdsk, read_eqdsk can be used'];
% data loaded to create eqdsks, can be used in geteqdskJET for other times to avoid to reload all
gdat_data.efitdata = efitdata;
-
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'mhd'}
params_eff = gdat_data.gdat_params;
@@ -699,7 +699,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
gdat_data.grids_1d = gdat_data2.grids_1d;
clear aa gdat_data2
end
- % load te if 'nete' was asked
+ % load te if 'nete' was asked
if ~isempty(data_request_effi{2})
params_eff.data_request = {'ppf',params_eff.source,data_request_effi{2}};
aa = gdat_jet(shot,params_eff);
@@ -878,7 +878,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
end
end
end
-
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ni','ti'}
nodenameeff_rho = 'rho';
@@ -939,7 +939,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
gdat_data.(data_request_eff).rhopol = gdat_data.x;
gdat_data.dimunits{1} = 'rhopol';
end
-
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ne_rho', 'te_rho', 'nete_rho'}
% if nete_rho, do first ne, then Te later (so fir stuff already done)
@@ -1215,7 +1215,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
% always start from ohmic so can use this time as base time since should yield full shot
fields_to_copy = {'data','units','dim','dimunits','t','x','data_fullpath','label','help','gdat_params', ...
- 'rad_bulk_h','rad_bulk_u','rad_bulk_avg'};
+ 'rad_bulk_h','rad_bulk_u','rad_bulk_avg','rad_bulk_t'};
fields_to_not_copy = {'shot','gdat_request'};
% total of each source in .data, but full data in subfield like pgyro in .ec, to check for nbi
params_eff = gdat_data.gdat_params;
@@ -1280,6 +1280,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
if any(strmatch('ic',gdat_data.gdat_params.source))
% ic
params_eff.data_request={'ppf','icrh','ptot'};
+ params_eff.data_request={'ppf','rff','ptot'};
try
ic=gdat_jet(shot,params_eff);
catch
@@ -1363,7 +1364,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
gdat_data.data = rad.data;
gdat_data.units = 'W';
gdat_data.data_fullpath = params_eff.data_request;
-
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'psi_edge'}
% psi at edge, 0 by construction in Liuqe, thus not given
@@ -1376,7 +1377,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end
return
- end
+ end
gdat_data.data = tracetdi.data.*0;
gdat_data.dim = tracetdi.dim;
gdat_data.t = gdat_data.dim{1};
@@ -1412,7 +1413,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
gdat_data.request_description = 'q(rhopol\_norm)';
% add grids_1d to have rhotor, etc
gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose);
-
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'rhotor', 'rhotor_edge', 'rhotor_norm', 'rhovol', 'volume_rho'}
% Phi(LCFS) = int(Bphi dSphi), can use Eq.(11) of "Tokamak coordinate conventions: COCOS" paper:
@@ -1437,7 +1438,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
phi_edge_2d = ones(length(x),1)*a(end,:);
gdat_data.rhotornorm = sqrt(phi./phi_edge_2d);
gdat_data.phi_tor = phi;
- params_eff = gdat_data.gdat_params;
+ params_eff = gdat_data.gdat_params;
params_eff.data_request='b0';
b0=gdat_jet(shot,params_eff);
gdat_data.b0 = interpos(b0.t,b0.data,t,-0.1);
@@ -1482,7 +1483,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
tracename2 = [nodenameeff2{2} '/' nodenameeff2{3}];
[a2,x2,t2,d2,e2]=rda_jet(shot,ppftype2,tracename2);
if e==0 || e2==0
- if gdat_params.nverbose>=3;
+ if gdat_params.nverbose>=3;
disp(['error after rda_jet in signal with nodenameeff= ' nodenameeff]);
disp(['or with nodenameeff2= ' nodenameeff2]);
disp(['cannot get ' data_request_eff]);
@@ -1511,7 +1512,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
return
end
-
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'rhovol','volume_rho','volume'}
% volume_rho = vol(rho); volume = vol(LCFS) = vol(rho=1);
@@ -1521,7 +1522,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
nodenameeff=[begstr 'vol' substr_liuqe];
end
tracetdi=tdi(nodenameeff);
- if isempty(tracetdi.data) || isempty(tracetdi.dim)
+ if isempty(tracetdi.data) || isempty(tracetdi.dim)
% try to run psitbxput
psitbxput_version = 1.3;
psitbxput(psitbxput_version,shot);
@@ -1558,7 +1559,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
end
end
gdat_data.t = gdat_data.dim{mapping_for_jet.gdat_timedim};
-
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'sxr'}
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
@@ -1656,7 +1657,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
else
tension = params_eff.tension;
end
-
+
aa=gdat_jet(shot,params_eff);
gdat_data.gdat_params = rmfield(params_eff,'data_request');
if ~isempty(params_eff.channel_index) && length(aa.data)~=numel(aa.data)
@@ -1677,7 +1678,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
error_status=901;
return
end
-
+
else
if (gdat_params.nverbose>=1); warning(['JET method=' mapping_for_jet.method ' not known yet, contact Olivier.Sauter@epfl.ch']); end
error_status=602;
diff --git a/crpptbx/JET/jet_requests_mapping.m b/crpptbx/JET/jet_requests_mapping.m
index 7e3cce87..ee25bb3a 100644
--- a/crpptbx/JET/jet_requests_mapping.m
+++ b/crpptbx/JET/jet_requests_mapping.m
@@ -48,6 +48,9 @@ switch lower(data_request)
mapping.timedim = 1;
mapping.method = 'signal';
mapping.expression = [{'ppf'},{'MAGN'},{'BVAC'}];
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request={''ppf'',''MAGN'',''BVAC''}; ' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);gdat_tmp.r0=2.96;'];
case 'beta'
mapping.label = 'beta\_tor diag [%]';
mapping.timedim = 1;
@@ -229,6 +232,26 @@ switch lower(data_request)
mapping.method = 'expression';
mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=[{''PPF''},{''EFIT''},{''FTOR''}];' ...
'gdat_tmp=gdat_jet(shot,params_eff);gdat_tmp.x=sqrt(gdat_tmp.x);gdat_tmp.dim{1}=gdat_tmp.x;'];
+ case {'p_lh', 'plh'}
+ mapping.label = 'p\_lh threshold';
+ mapping.timedim = 1;
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''nel'';' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);params_eff.data_request=''b0'';gdat_tmpb0=gdat_jet(shot,params_eff);' ...
+ 'params_eff.data_request=''surface_edge'';gdat_tmp_s=gdat_jet(shot,params_eff);' ...
+ 's_nel=interpos(gdat_tmp_s.t,gdat_tmp_s.data,gdat_tmp.t,-1.);b0_nel=interpos(gdat_tmpb0.t,gdat_tmpb0.data,gdat_tmp.t,-1.);' ...
+ 'gdat_tmp.data = 0.0488e6.*(gdat_tmp.data/1e20).^0.717.*abs(b0_nel).^0.803.*s_nel.^0.941;' ...
+ 'gdat_tmp.label=''P\_LH threshold'';'];
+ case {'p_lh_a_r', 'plh_a_r'}
+ mapping.label = 'p\_lh threshold';
+ mapping.timedim = 1;
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''nel'';' ...
+ 'gdat_tmp=gdat_jet(shot,params_eff);params_eff.data_request=''b0'';gdat_tmpb0=gdat_jet(shot,params_eff);' ...
+ 'params_eff.data_request=''a_minor'';gdat_tmp_a=gdat_jet(shot,params_eff);' ...
+ 'a_nel=interpos(gdat_tmp_a.t,gdat_tmp_a.data,gdat_tmp.t,-1.);b0_nel=interpos(gdat_tmpb0.t,gdat_tmpb0.data,gdat_tmp.t,-1.);' ...
+ 'gdat_tmp.data = 2.15e6.*abs(gdat_tmp.data/1e20).^0.782.*abs(b0_nel).^0.772.*abs(a_nel).^0.975.*abs(gdat_tmpb0.r0).^0.999;' ...
+ 'gdat_tmp.label=''P\_LH threshold'';'];
case {'p_ohmic', 'p_ohm', 'pohm'}
mapping.label = 'p\_ohmic';
mapping.timedim = 1;
@@ -377,7 +400,7 @@ switch lower(data_request)
mapping.label = data_request;
mapping.method = 'signal'; % assume a full tracename is given, so just try with tdi (could check there are some ":", etc...)
mapping.expression = data_request;
-
+
end
if isempty(mapping.gdat_timedim)
diff --git a/crpptbx/TCV/tcv_requests_mapping.m b/crpptbx/TCV/tcv_requests_mapping.m
index 5195dd47..1dd71bbe 100644
--- a/crpptbx/TCV/tcv_requests_mapping.m
+++ b/crpptbx/TCV/tcv_requests_mapping.m
@@ -255,6 +255,26 @@ switch lower(data_request)
% node not filled in and in any case need special case for
mapping.label = 'toroidal_flux';
mapping.method = 'switchcase';
+ case {'p_lh', 'plh'}
+ mapping.label = 'p\_lh threshold';
+ mapping.timedim = 1;
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''nel'';' ...
+ 'gdat_tmp=gdat_tcv(shot,params_eff);params_eff.data_request=''b0'';gdat_tmpb0=gdat_tcv(shot,params_eff);' ...
+ 'params_eff.data_request=''surface_edge'';gdat_tmp_s=gdat_tcv(shot,params_eff);' ...
+ 's_nel=interpos(gdat_tmp_s.t,gdat_tmp_s.data,gdat_tmp.t,-1.);b0_nel=interpos(gdat_tmpb0.t,gdat_tmpb0.data,gdat_tmp.t,-1.);' ...
+ 'gdat_tmp.data = 0.0488e6.*(gdat_tmp.data/1e20).^0.717.*abs(b0_nel).^0.803.*s_nel.^0.941;' ...
+ 'gdat_tmp.label=''P\_LH threshold'';'];
+ case {'p_lh_a_r', 'plh_a_r'}
+ mapping.label = 'p\_lh threshold';
+ mapping.timedim = 1;
+ mapping.method = 'expression';
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''nel'';' ...
+ 'gdat_tmp=gdat_tcv(shot,params_eff);params_eff.data_request=''b0'';gdat_tmpb0=gdat_tcv(shot,params_eff);' ...
+ 'params_eff.data_request=''a_minor'';gdat_tmp_a=gdat_tcv(shot,params_eff);' ...
+ 'a_nel=interpos(gdat_tmp_a.t,gdat_tmp_a.data,gdat_tmp.t,-1.);b0_nel=interpos(gdat_tmpb0.t,gdat_tmpb0.data,gdat_tmp.t,-1.);' ...
+ 'gdat_tmp.data = 2.15e6.*(gdat_tmp.data/1e20).^0.782.*abs(b0_nel).^0.772.*a_nel.^0.975.*gdat_tmpb0.r0.^0.999;' ...
+ 'gdat_tmp.label=''P\_LH threshold'';'];
case 'psi'
mapping.timedim = 1;
mapping.label = 'psi(R,Z)';
@@ -369,6 +389,17 @@ switch lower(data_request)
mapping.label = 'rtc\_signals';
mapping.method = 'switchcase';
mapping.expression = '';
+ case 'surface_rho'
+ mapping.timedim = 1;
+ mapping.label = 'surface(rho)';
+ mapping.method = 'tdiliuqe';
+ mapping.expression = 'tcv_eq(''''surf'''',''''LIUQE.M'''')';
+ 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'''''''')'';' ...
+ 'gdat_tmp=gdat_tcv(shot,params_eff); gdat_tmp.dim = {gdat_tmp.t}; gdat_tmp.x=[]; gdat_tmp.data= gdat_tmp.data(end,:);'];
case 'sxr'
mapping.timedim = 1;
mapping.gdat_timedim = 2;
--
GitLab