From f8cdc4986ae6052f3e935efc2aa2d55dcfbfe8b4 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Wed, 5 Jul 2017 08:39:55 +0000
Subject: [PATCH] added powers to jet and radiated from bolo
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@7760 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/JET/gdat_jet.m | 364 +++++++++++++++++++++++------
crpptbx/JET/jet_requests_mapping.m | 9 +-
crpptbx/gdat_plot.m | 2 +-
3 files changed, 294 insertions(+), 81 deletions(-)
diff --git a/crpptbx/JET/gdat_jet.m b/crpptbx/JET/gdat_jet.m
index d394cce2..925ec592 100644
--- a/crpptbx/JET/gdat_jet.m
+++ b/crpptbx/JET/gdat_jet.m
@@ -782,12 +782,123 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ne','te'}
+ nodenameeff_rho = 'rho';
+ params_eff = gdat_data.gdat_params;
+ if isfield(params_eff,'source') && ~isempty(params_eff.source)
+ if strcmp(lower(params_eff.source),'chain2')
+ params_eff.source = 'hrt2';
+ end
+ if strcmp(lower(params_eff.source),'hrt2')
+ nodenameeff_rho = []; % profiles already on rho
+ end
+ else
+ params_eff.source = 'hrts';
+ end
+ gdat_data.gdat_params = params_eff;
+ params_eff.doplot = 0;
% ne or Te from Thomson data on raw z mesh vs (z,t)
- 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;
+ params_eff.data_request = {'ppf',params_eff.source,data_request_eff};
+ aa = gdat_jet(shot,params_eff);
+ if isempty(aa.data) || isempty(aa.t)
+ if gdat_params.nverbose>=3;
+ aa
+ disp(['with data_request= ' data_request_eff])
+ end
+ return
+ end
+ gdat_data.data = aa.data;
+ gdat_data.dim = aa.dim;
+ gdat_data.t = aa.dim{2};
+ gdat_data.x = aa.dim{1};
+ gdat_data.dimunits=[{'R [m]'} ; {'time [s]'}];
+ gdat_data.data_fullpath=[data_request_eff ' from ' params_eff.source];
+ gdat_data.(data_request_eff).data = aa.data;
+ gdat_data.(data_request_eff).t = aa.t;
+ gdat_data.(data_request_eff).r = aa.x;
+ if any(strcmp(fieldnames(aa),'units')) && ~isempty(aa.units)
+ gdat_data.units=aa.units;
+ else
+ if strcmp(data_request_eff,'ne')
+ gdat_data.units = 'm^{-3}';
+ elseif strcmp(data_request_eff,'te')
+ gdat_data.units = 'eV';
+ end
+ end
+ params_eff.data_request = {'ppf',params_eff.source,['d' data_request_eff]};
+ aaerr = gdat_jet(shot,params_eff);
+ gdat_data.error_bar = aaerr.data;
+ gdat_data.(data_request_eff).error_bar = aaerr.data;
+ if ~isempty(nodenameeff_rho)
+ params_eff.data_request = {'ppf',params_eff.source,'rho'};
+ aarho = gdat_jet(shot,params_eff);
+ gdat_data.rhopol = abs(aarho.data);
+ gdat_data.(data_request_eff).rhopol = abs(aarho.data); % keep rho >0
+ else
+ gdat_data.rhopol = gdat_data.x;
+ gdat_data.(data_request_eff).rhopol = gdat_data.x;
+ gdat_data.dimunits{1} = 'rhopol';
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'ni','ti'}
+ nodenameeff_rho = 'rho';
+ data_type_eff = data_request_eff;
+ if strcmp(data_request_eff,'ni'); data_type_eff = 'dd'; end
+ params_eff = gdat_data.gdat_params;
+ if isfield(params_eff,'source') && ~isempty(params_eff.source)
+ if strcmp(lower(params_eff.source),'chain2')
+ params_eff.source = [data_request_eff(1) 'ion'];
+ end
+ if strcmp(lower(params_eff.source(2:end)),'ion')
+ nodenameeff_rho = []; % profiles already on rho
+ end
+ else
+ params_eff.source = [datar_equest_eff(1) 'ion']; % only chain2 at this stage
+ end
+ gdat_data.gdat_params = params_eff;
+ params_eff.doplot = 0;
+ % ne or Te from Thomson data on raw z mesh vs (z,t)
+ params_eff.data_request = {'ppf',params_eff.source,data_type_eff};
+ aa = gdat_jet(shot,params_eff);
+ if isempty(aa.data) || isempty(aa.t)
+ if gdat_params.nverbose>=3;
+ aa
+ disp(['with data_request= ' data_request_eff])
+ end
+ return
+ end
+ gdat_data.data = aa.data;
+ gdat_data.dim = aa.dim;
+ gdat_data.t = aa.dim{2};
+ gdat_data.x = aa.dim{1};
+ gdat_data.dimunits=[{'R [m]'} ; {'time [s]'}];
+ gdat_data.data_fullpath=[data_request_eff ' from ' params_eff.source];
+ gdat_data.(data_request_eff).data = aa.data;
+ gdat_data.(data_request_eff).t = aa.t;
+ gdat_data.(data_request_eff).r = aa.x;
+ if any(strcmp(fieldnames(aa),'units')) && ~isempty(aa.units)
+ gdat_data.units=aa.units;
+ else
+ if strcmp(data_request_eff,'ni')
+ gdat_data.units = 'm^{-3}';
+ elseif strcmp(data_request_eff,'ti')
+ gdat_data.units = 'eV';
+ end
+ end
+ params_eff.data_request = {'ppf',params_eff.source,['e' data_type_eff]}; % not given for dd but ok will be empty
+ aaerr = gdat_jet(shot,params_eff);
+ gdat_data.error_bar = aaerr.data;
+ gdat_data.(data_request_eff).error_bar = aaerr.data;
+ if ~isempty(nodenameeff_rho)
+ params_eff.data_request = {'ppf',params_eff.source,'rho'};
+ aarho = gdat_jet(shot,params_eff);
+ gdat_data.rhopol = abs(aarho.data);
+ gdat_data.(data_request_eff).rhopol = abs(aarho.data); % keep rho >0
+ else
+ gdat_data.rhopol = gdat_data.x;
+ gdat_data.(data_request_eff).rhopol = gdat_data.x;
+ gdat_data.dimunits{1} = 'rhopol';
end
- [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,gdat_data.gdat_params.edge,gdat_params.nverbose);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ne_rho', 'te_rho', 'nete_rho'}
@@ -1023,83 +1134,151 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
case {'powers'}
% note: same time array for all main, ec, ohm, nbi, ...
% At this stage fill just ech, later add nbi
-
- % EC
- nodenameeff='\results::toray.input:p_gyro';
- tracetdi=tdi(nodenameeff);
- if isempty(tracetdi.data) || isempty(tracetdi.dim) % || ischar(tracetdi.data) (to add?)
- if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
- gdat_data.ec.data = [];
- gdat_data.ec.units = [];
- gdat_data.ec.dim=[];
- gdat_data.ec.dimunits=[];
- gdat_data.ec.t=[];
- gdat_data.ec.x=[];
- gdat_data.ec.data_fullpath=[];
- gdat_data.ec.label='';
+ sources_avail = {'ohm','nbi','ic','rad'}; % note should allow ech, nbh, ohmic in parameter sources
+ for i=1:length(sources_avail)
+ gdat_data.(sources_avail{i}).data = [];
+ gdat_data.(sources_avail{i}).units = [];
+ gdat_data.(sources_avail{i}).dim=[];
+ gdat_data.(sources_avail{i}).dimunits=[];
+ gdat_data.(sources_avail{i}).t=[];
+ gdat_data.(sources_avail{i}).x=[];
+ gdat_data.(sources_avail{i}).data_fullpath=[];
+ gdat_data.(sources_avail{i}).label=[];
+ end
+ if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
+ gdat_data.gdat_params.source = sources_avail;
+ elseif ~iscell(gdat_data.gdat_params.source)
+ if ischar(gdat_data.gdat_params.source)
+ gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source);
+ if ~any(strmatch(gdat_data.gdat_params.source,lower(sources_avail)))
+ if (gdat_params.nverbose>=1)
+ warning(['source= ' gdat_data.gdat_params.source ' is not part of the available sources: ' sprintf('''%s'' ',sources_avail{:})]);
+ end
+ return
+ else
+ gdat_data.gdat_params.source = {gdat_data.gdat_params.source};
+ end
+ else
+ if (gdat_params.nverbose>=1); warning([' source parameter not compatible with: ' sprintf('''%s'' ',sources_avail{:})]); end
+ return
+ end
else
- 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';
+ for i=1:length(gdat_data.gdat_params.source)
+ gdat_data.gdat_params.source{i} = lower(gdat_data.gdat_params.source{i});
+ if ~any(strmatch(gdat_data.gdat_params.source{i},lower(sources_avail)))
+ if gdat_data.gdat_params.nverbose>=1
+ warning(['source = ' gdat_data.gdat_params.source{i} ' not expected with data_request= ' data_request_eff])
+ end
+ end
+ end
end
+ % always start from ohmic so can use this time as base time since should yield full shot
- % Ohmic
- gdat_data.ohm.data = [];
- gdat_data.ohm.units = '';
- gdat_data.ohm.dim=gdat_data.dim;
- gdat_data.ohm.dimunits=gdat_data.dimunits;
- gdat_data.ohm.t=[];
- gdat_data.ohm.x=[];
- gdat_data.ohm.data_fullpath=[];
- gdat_data.ohm.label='';
- % get ohmic power simply from vloop*Ip (minus sign for JET)
- ip=gdat([],'ip');
- vloop=gdat([],'vloop');
- tension = -1e5;
- if isempty(gdat_data.t)
- gdat_data.t = vloop.t;
- gdat_data.dim{1} = gdat_data.t;
- gdat_data.ohm.data = vloop.data;
- gdat_data.dimunits{1} = 's';
- gdat_data.units = 'W';
- else
- 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;
- end
- gdat_data.ohm.units = gdat_data.units;
- 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 each and total
- gdat_data.data(:,1) = gdat_data.ohm.data;
- if ~isempty(gdat_data.ec.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_{ohm};P_{EC};P_{tot}';
+ fields_to_copy = {'data','units','dim','dimunits','t','x','data_fullpath','label','help','gdat_params'};
+ 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;
+ params_eff.source=[]; % use default for individual calls
+
+ % ohmic, use its time-base
+ params_eff.data_request='p_ohmic';
+ try
+ ohm=gdat_jet(shot,params_eff);
+ catch
+ ohm.data = [];
+ ohm.dim = [];
+ end
+ if ~isempty(ohm.data) && ~isempty(ohm.dim)
+ for i=1:length(fields_to_copy)
+ if isfield(ohm,fields_to_copy{i})
+ gdat_data.ohm.(fields_to_copy{i}) = ohm.(fields_to_copy{i});
+ end
+ end
+ gdat_data.ohm.raw_data = gdat_data.ohm.data;
else
- gdat_data.dim{2} = [1];
- gdat_data.dimunits{2} = 'Pohm=Ptot';
- gdat_data.data_fullpath=['tot power from ohm'];
- gdat_data.label = 'P_{tot}=P_{ohm}';
+ if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
+ return
end
+ mapping_for_jet.timedim = 1; mapping_for_jet.gdat_timedim = 1;
+ taus = -10;
+ %
+ % add each source in main.data, on ohm time array
+ gdat_data.t = gdat_data.ohm.t;
+ gdat_data.dim{1} = gdat_data.t;
+ gdat_data.dimunits{1} = 's';
+ gdat_data.ohm.data = interpos(gdat_data.t,gdat_data.ohm.raw_data,5.*taus);
+ gdat_data.data = reshape(gdat_data.ohm.data,length(gdat_data.t),1);
+ gdat_data.ohm.tension = 5.*taus;
+ gdat_data.x =[1];
+ gdat_data.label={'P_{ohm}'};
+ gdat_data.units = 'W';
+ %
+ if any(strmatch('nb',gdat_data.gdat_params.source))
+ % nbi
+ params_eff.data_request={'ppf','nbi','ptot'};
+ try
+ nbi=gdat_jet(shot,params_eff);
+ catch
+ end
+ if ~isempty(nbi.data) && ~isempty(nbi.dim)
+ for i=1:length(fields_to_copy)
+ if isfield(nbi,fields_to_copy{i})
+ gdat_data.nbi.(fields_to_copy{i}) = nbi.(fields_to_copy{i});
+ end
+ end
+ % add to main with linear interpolation and 0 for extrapolated values
+ gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi.t,gdat_data.nbi.data(:,end),gdat_data.t);
+ gdat_data.x(end+1) =gdat_data.x(end)+1;
+ gdat_data.label{end+1}='P_{nbi}';
+ end
+ end
+ %
+ if any(strmatch('ic',gdat_data.gdat_params.source))
+ % ic
+ params_eff.data_request={'ppf','icrh','ptot'};
+ try
+ ic=gdat_jet(shot,params_eff);
+ catch
+ end
+ if ~isempty(ic.data) && ~isempty(ic.dim)
+ for i=1:length(fields_to_copy)
+ if isfield(ic,fields_to_copy{i})
+ gdat_data.ic.(fields_to_copy{i}) = ic.(fields_to_copy{i});
+ end
+ end
+ gdat_data.data(:,end+1) = interpos(-21,gdat_data.ic.t,gdat_data.ic.data,gdat_data.t);
+ gdat_data.x(end+1) =gdat_data.x(end)+1;
+ gdat_data.label{end+1}='P_{ic}';
+ end
+ end
+ if any(strmatch('rad',gdat_data.gdat_params.source))
+ % radiated power
+ params_eff.data_request='p_rad';
+ try
+ rad=gdat_jet(shot,params_eff);
+ catch
+ end
+ if ~isempty(rad.data) && ~isempty(rad.dim)
+ for i=1:length(fields_to_copy)
+ if isfield(rad,fields_to_copy{i})
+ gdat_data.rad.(fields_to_copy{i}) = rad.(fields_to_copy{i});
+ end
+ end
+ % add to main with linear interpolation and 0 for extrapolated values
+ gdat_data.data(:,end+1) = interpos(-21,gdat_data.rad.t,gdat_data.rad.data,gdat_data.t);
+ gdat_data.x(end+1) =gdat_data.x(end)+1;
+ gdat_data.label{end+1}='P_{rad}';
+ end
+ end
+ % add tot power
+ sources_coeff = [ 1, 1, 1, 0];
+ aa=sum(gdat_data.data.*repmat(reshape(sources_coeff,1,numel(sources_coeff)),size(gdat_data.data,1),1),2);
+ gdat_data.data(:,end+1) = aa;
+ % gdat_data.data(:,end+1) = sum(gdat_data.data,2);
+ gdat_data.label{end+1}='P_{tot} without Prad';
+ gdat_data.x(end+1) =gdat_data.x(end)+1;
+ gdat_data.dim{2} = gdat_data.x; gdat_data.dimunits{2} = '';
+ gdat_data.data_fullpath = 'tot powers from each sources, and total power in .data(:,end)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'q_rho'}
@@ -1128,6 +1307,39 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
% add grids_1d to have rhotor, etc
gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose);
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'p_rad', 'prad'}
+ params_eff = gdat_data.gdat_params;
+ if ~isfield(params_eff,'source') || isempty(params_eff.source)
+ params_eff.source = 'bolo';
+ end
+ gdat_data.gdat_params = params_eff;
+ switch params_eff.source
+ case 'prad'
+ params_eff.data_request = {'ppf','prad','prad'};
+ rad=gdat_jet(shot,params_eff);
+ gdat_data.label={'P_{rad} from prad'};
+ otherwise
+ % 'bolo' by default
+ gdat_data.gdat_params.source = 'bolo';
+ params_eff.data_request = {'ppf','bolo','topi'};
+ rad=gdat_jet(shot,params_eff);
+ gdat_data.label={'P_{rad} from bolo/topi'};
+ params_eff.data_request = {'ppf','bolo','tobh'};
+ rad_bulk1=gdat_jet(shot,params_eff);
+ gdat_data.rad_bulk_h = rad_bulk1.data;
+ params_eff.data_request = {'ppf','bolo','tobu'};
+ rad_bulk2=gdat_jet(shot,params_eff);
+ gdat_data.rad_bulk_u = rad_bulk2.data;
+ gdat_data.rad_bulk_avg = 0.5.*(rad_bulk1.data+rad_bulk2.data);
+ end
+ gdat_data.t = rad.t;
+ gdat_data.dim{1} = gdat_data.t;
+ gdat_data.dimunits{1} = 's';
+ 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
diff --git a/crpptbx/JET/jet_requests_mapping.m b/crpptbx/JET/jet_requests_mapping.m
index f843d20f..03bde4ee 100644
--- a/crpptbx/JET/jet_requests_mapping.m
+++ b/crpptbx/JET/jet_requests_mapping.m
@@ -129,7 +129,7 @@ switch lower(data_request)
case 'mhd'
mapping.timedim = 1;
mapping.label = 'n=1 and n=2 amplitude';
- mapping.method = 'expression';
+ mapping.method = 'expression'; % {'JPF','DA','C1M-H306'} {'JPF','DA','C1M-T009'} (big) {'JPF','DA','C1M-MIT27'}
mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request={''JPF'',''DA'',''C1H-I301''}; ' ...
'gdat_tmp=gdat_jet(shot,params_eff);gdat_tmp.data=reshape(gdat_tmp.data,length(gdat_tmp.data),1 );' ...
'gdat_tmp.dim{1}=gdat_tmp.t;gdat_tmp.dim{2}=[1 2];gdat_tmp.x=gdat_tmp.dim{2};' ...
@@ -203,7 +203,7 @@ switch lower(data_request)
mapping.timedim = 1;
mapping.method = 'signal';
mapping.expression = [{'PPF'},{'EQUI'},{'PHIT'}]; % note this is only chain2, so should check with efit...
- case 'p_ohmic'
+ case {'p_ohmic', 'p_ohm', 'pohm'}
mapping.label = 'p\_ohmic';
mapping.timedim = 1;
mapping.method = 'signal';
@@ -211,8 +211,9 @@ switch lower(data_request)
case {'p_rad', 'prad'}
mapping.label = 'p\_rad';
mapping.timedim = 1;
- mapping.method = 'signal';
- mapping.expression = [{'PPF'},{'PRAD'},{'PRAD'}];
+% $$$ mapping.method = 'signal';
+% $$$ mapping.expression = [{'PPF'},{'bolo'},{'topi'}];
+ mapping.method = 'switchcase';
case 'powers'
mapping.timedim = 1;
mapping.label = 'various powers';
diff --git a/crpptbx/gdat_plot.m b/crpptbx/gdat_plot.m
index 5cdee4d7..3377affc 100644
--- a/crpptbx/gdat_plot.m
+++ b/crpptbx/gdat_plot.m
@@ -91,7 +91,7 @@ if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(
if strcmp(gdat_data.gdat_request,'mhd')
% special case, add legend instead of ylabel
delete(hylabel);
- legend(gdat_data.label,2);
+ legend(gdat_data.label);
% add spectrogram per signal
mhd_sum_data = 0.;
nfft=1024;
--
GitLab