From ea04520378e533de1631b8a6decdf3fad101c047 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Fri, 7 Jul 2017 13:42:03 +0000
Subject: [PATCH] added ne, te and nete including fit with edge value, minimum
value and mask options, latter required since chain2 does not do it
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@7820 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/JET/gdat_jet.m | 241 +++++++++++++++++++----------
crpptbx/JET/jet_help_parameters.m | 1 +
crpptbx/JET/jet_requests_mapping.m | 4 +
3 files changed, 164 insertions(+), 82 deletions(-)
diff --git a/crpptbx/JET/gdat_jet.m b/crpptbx/JET/gdat_jet.m
index 68b3ff4d..3dfdbe76 100644
--- a/crpptbx/JET/gdat_jet.m
+++ b/crpptbx/JET/gdat_jet.m
@@ -52,6 +52,9 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_jet(shot,data_req
% For JET, the specific trace can be given as:
% aa==gdat(51994,[{'ppf'},{'efit'},{'xip'}],'machine','jet','doplot',1);
%
+% Example to mask some hrts channels in the fit:
+% nete=gdat(92442,'nete','machine','jet','doplot',1,'fit_mask',[25:27]);
+%
% Comments for local developer:
% This gdat is just a "header routine" calling the gdat for the specific machine gdat_`machine`.m which can be called
% directly, thus which should be able to treat the same type of input arguments
@@ -613,7 +616,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- case {'ne','te'}
+ case {'ne','te', 'nete'}
nodenameeff_rho = 'rho';
if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit) || ~isnumeric(gdat_data.gdat_params.fit)
gdat_data.gdat_params.fit = 1; % default do fit
@@ -637,8 +640,15 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
end
params_eff = gdat_data.gdat_params;
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_request_eff};
+ % ne and/or Te from Thomson data on raw z mesh vs (z,t)
+ data_request_effi{1} = data_request_eff;
+ if strcmp(data_request_eff,'nete')
+ data_request_effi{1} = 'ne'; % start with ne
+ data_request_effi{2} = 'te';
+ else
+ data_request_effi{2} = [];
+ end
+ params_eff.data_request = {'ppf',params_eff.source,data_request_effi{1}};
aa = gdat_jet(shot,params_eff);
if isempty(aa.data) || isempty(aa.t)
if gdat_params.nverbose>=3;
@@ -653,30 +663,31 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
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;
+ gdat_data.(data_request_effi{1}).data = aa.data;
+ gdat_data.(data_request_effi{1}).t = aa.t;
+ gdat_data.(data_request_effi{1}).r = aa.x;
if any(strcmp(fieldnames(aa),'units')) && ~isempty(aa.units)
gdat_data.units=aa.units;
else
- if strcmp(data_request_eff,'ne')
+ if strcmp(data_request_effi{1},'ne')
gdat_data.units = 'm^{-3}';
- elseif strcmp(data_request_eff,'te')
+ elseif strcmp(data_request_effi{1},'te')
gdat_data.units = 'eV';
end
end
- params_eff.data_request = {'ppf',params_eff.source,['d' data_request_eff]};
+ params_eff.data_request = {'ppf',params_eff.source,['d' data_request_effi{1}]};
aaerr = gdat_jet(shot,params_eff);
gdat_data.error_bar = aaerr.data;
- gdat_data.(data_request_eff).error_bar = aaerr.data;
+ gdat_data.(data_request_effi{1}).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
+ gdat_data.(data_request_effi{1}).rhopol = abs(aarho.data); % keep rho >0
+ gdat_data.(data_request_effi{1}).rhopol_sign = aarho.data; % to be able to distinguish channels
else
gdat_data.rhopol = gdat_data.x;
- gdat_data.(data_request_eff).rhopol = gdat_data.x;
+ gdat_data.(data_request_effi{1}).rhopol = gdat_data.x;
gdat_data.dimunits{1} = 'rhopol';
end
if length(gdat_data.x) == numel(gdat_data.rhopol)
@@ -688,18 +699,49 @@ 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
+ if ~isempty(data_request_effi{2})
+ params_eff.data_request = {'ppf',params_eff.source,data_request_effi{2}};
+ 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_effi{2}])
+ end
+ return
+ end
+ gdat_data.(data_request_effi{2}).data = aa.data;
+ gdat_data.(data_request_effi{2}).t = aa.t;
+ gdat_data.(data_request_effi{2}).r = aa.x;
+ if any(strcmp(fieldnames(aa),'units')) && ~isempty(aa.units)
+ gdat_data.(data_request_effi{2}).units=aa.units;
+ else
+ if strcmp(data_request_effi{2},'ne')
+ gdat_data.(data_request_effi{2}).units = 'm^{-3}';
+ elseif strcmp(data_request_effi{2},'te')
+ gdat_data.(data_request_effi{2}).units = 'eV';
+ end
+ end
+ params_eff.data_request = {'ppf',params_eff.source,['d' data_request_effi{2}]};
+ aaerr = gdat_jet(shot,params_eff);
+ gdat_data.(data_request_effi{2}).error_bar = aaerr.data;
+ gdat_data.(data_request_effi{2}).rhopol = gdat_data.(data_request_effi{1}).rhopol;
+ gdat_data.(data_request_effi{2}).rhopol_sign = gdat_data.(data_request_effi{1}).rhopol_sign;
+ % construct pressure in .data
+ gdat_data.data = 1.602e-19.* gdat_data.(data_request_effi{1}).data .* gdat_data.(data_request_effi{2}).data;
+ end
% defaults for fits, so user always gets std structure
gdat_data.fit.rhotornorm = []; % same for both ne and te
gdat_data.fit.rhopolnorm = [];
gdat_data.fit.t = [];
gdat_data.fit.te.data = [];
- gdat_data.fit.te.drhotornorm = [];
+ gdat_data.fit.te.d_over_drhotornorm = [];
gdat_data.fit.ne.data = [];
- gdat_data.fit.ne.drhotornorm = [];
+ gdat_data.fit.ne.d_over_drhotornorm = [];
gdat_data.fit.raw.te.data = [];
- gdat_data.fit.raw.te.rhotornorm = [];
+ gdat_data.fit.raw.te.rhopolnorm = [];
gdat_data.fit.raw.ne.data = [];
- gdat_data.fit.raw.ne.rhotornorm = [];
+ gdat_data.fit.raw.ne.rhopolnorm = [];
fit_tension_default = -1;
if isfield(gdat_data.gdat_params,'fit_tension')
fit_tension = gdat_data.gdat_params.fit_tension;
@@ -721,6 +763,27 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
fit_nb_rho_points = 201;
end
gdat_data.gdat_params.fit_nb_rho_points = fit_nb_rho_points;
+ if isfield(gdat_data.gdat_params,'fit_ne_min')
+ fit_ne_min = gdat_data.gdat_params.fit_ne_min;
+ else
+ fit_ne_min = 1e17;
+ end
+ if isfield(gdat_data.gdat_params,'fit_te_min')
+ fit_te_min = gdat_data.gdat_params.fit_te_min;
+ else
+ fit_te_min = 5;
+ end
+ gdat_data.gdat_params.fit_te_min = fit_te_min;
+ fit_mask = -1; % we mask on the indices, thus -1 means no mask
+ if isfield(gdat_data.gdat_params,'fit_mask') && isnumeric(gdat_data.gdat_params.fit_mask) ...
+ && ~isempty(gdat_data.gdat_params.fit_mask)
+ % channel index to mask
+ fit_mask = gdat_data.gdat_params.fit_mask;
+ else
+ fit_mask = fit_mask;
+ end
+ gdat_data.gdat_params.fit_mask = fit_mask;
+ gdat_data.fit.mask = fit_mask;
%
if gdat_data.gdat_params.fit==1
% add fits
@@ -728,77 +791,91 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
rhopolfit = linspace(0,1,fit_nb_rho_points)';
gdat_data.fit.rhotornorm = NaN*ones(length(rhopolfit),length(gdat_data.fit.t));
gdat_data.fit.rhopolnorm = rhopolfit;
- if any(strfind(data_request_eff,'te'))
- gdat_data.fit.raw.te.data = NaN*ones(size(gdat_data.te.data));
- gdat_data.fit.raw.te.error_bar = NaN*ones(size(gdat_data.te.data));
- gdat_data.fit.raw.te.rhopolnorm = NaN*ones(size(gdat_data.te.data));
- gdat_data.fit.te.data = gdat_data.fit.rhotornorm;
- gdat_data.fit.te.drhotornorm = gdat_data.fit.rhotornorm;
- end
- if any(strfind(data_request_eff,'ne'))
- gdat_data.fit.raw.ne.data = NaN*ones(size(gdat_data.ne.data));
- gdat_data.fit.raw.ne.error_bar = NaN*ones(size(gdat_data.ne.data));
- gdat_data.fit.raw.ne.rhopolnorm = NaN*ones(size(gdat_data.ne.data));
- gdat_data.fit.ne.data =gdat_data.fit.rhotornorm;
- gdat_data.fit.ne.drhotornorm = gdat_data.fit.rhotornorm;
- end
+ % common part between ne and te
+ indices_ok=[1:size(gdat_data.data,1)]';
+ indices_ok=setdiff(indices_ok,fit_mask);
for it=1:length(gdat_data.t)
- % make rhopol->rhotor transformation for each time since equilibrium might have changed
- irhook=find(gdat_data.grids_1d.rhopolnorm(:,it)>0 & gdat_data.grids_1d.rhopolnorm(:,it)<1); % no need for ~isnan
- [rhoeff isort]=sort(gdat_data.grids_1d.rhopolnorm(irhook,it));
- gdat_data.fit.rhotornorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.grids_1d.rhotornorm(irhook(isort),it); 1],rhopolfit,-0.1,[2 2],[0 1]);
- if any(strfind(data_request_eff,'te'))
- idatate = find(gdat_data.te.data(:,it)>1 & gdat_data.grids_1d.rhopolnorm(:,it)<=1.05);
- if length(idatate)>0
- gdat_data.fit.raw.te.rhopolnorm(idatate,it) = gdat_data.grids_1d.rhopolnorm(idatate,it);
- gdat_data.fit.raw.te.data(idatate,it) = gdat_data.te.data(idatate,it);
- gdat_data.fit.raw.te.error_bar(idatate,it) = gdat_data.te.error_bar(idatate,it);
- [rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhopolnorm(idatate,it));
- rhoeff = [0; rhoeff];
- teeff = gdat_data.te.data(idatate(irhoeff),it);
- te_err_eff = gdat_data.te.error_bar(idatate(irhoeff),it);
- % they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
- ij=find(teeff./te_err_eff>10.);
- if ~isempty(ij); te_err_eff(ij) = teeff(ij)./0.1; end
- %
- teeff = [teeff(1); teeff];
- te_err_eff = [1e4; te_err_eff];
- if gdat_data.gdat_params.fit_te_edge < 0
- [gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhopolnorm(:,it)] = interpos(rhoeff,teeff,rhopolfit,fit_tension.te,[1 0],[0 ...
- 0],te_err_eff);
- else
- [gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhopolnorm(:,it)] = ...
- interpos(rhoeff,teeff,rhopolfit,fit_tension.te,[1 2],[0 gdat_data.gdat_params.fit_te_edge],te_err_eff);
+ % make rhopol->rhotor transformation for each time since equilibrium might have changed
+ irhook=find(gdat_data.grids_1d.rhopolnorm(indices_ok,it)>0 & gdat_data.grids_1d.rhopolnorm(indices_ok,it)<1); % no need for ~isnan
+ [rhoeff isort]=sort(gdat_data.grids_1d.rhopolnorm(indices_ok(irhook),it));
+ gdat_data.fit.rhotornorm(:,it)=interpos([0; rhoeff; 1], ...
+ [0; gdat_data.grids_1d.rhotornorm(indices_ok(irhook(isort)),it); 1],rhopolfit,-0.1,[2 2],[0 1]);
+ end
+ for i=1:2
+ if strcmp(data_request_effi{i},'te')
+ gdat_data.fit.raw.te.data = NaN*ones(size(gdat_data.te.data));
+ gdat_data.fit.raw.te.error_bar = NaN*ones(size(gdat_data.te.data));
+ gdat_data.fit.raw.te.rhopolnorm = NaN*ones(size(gdat_data.te.data));
+ gdat_data.fit.te.data = gdat_data.fit.rhotornorm;
+ gdat_data.fit.te.d_over_drhotornorm = gdat_data.fit.rhotornorm;
+ elseif strcmp(data_request_effi{i},'ne')
+ gdat_data.fit.raw.ne.data = NaN*ones(size(gdat_data.ne.data));
+ gdat_data.fit.raw.ne.error_bar = NaN*ones(size(gdat_data.ne.data));
+ gdat_data.fit.raw.ne.rhopolnorm = NaN*ones(size(gdat_data.ne.data));
+ gdat_data.fit.ne.data =gdat_data.fit.rhotornorm;
+ gdat_data.fit.ne.d_over_drhotornorm = gdat_data.fit.rhotornorm;
+ end
+ for it=1:length(gdat_data.t)
+ if strcmp(data_request_effi{i},'te')
+ idatate = find(gdat_data.te.data(indices_ok,it)>1 & gdat_data.grids_1d.rhopolnorm(indices_ok,it)<=1.05);
+ idatate = indices_ok(idatate);
+ if length(idatate)>0
+ gdat_data.fit.raw.te.rhopolnorm(idatate,it) = gdat_data.grids_1d.rhopolnorm(idatate,it);
+ gdat_data.fit.raw.te.data(idatate,it) = gdat_data.te.data(idatate,it);
+ gdat_data.fit.raw.te.error_bar(idatate,it) = gdat_data.te.error_bar(idatate,it);
+ [rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhopolnorm(idatate,it));
+ rhoeff = [0; rhoeff];
+ teeff = gdat_data.te.data(idatate(irhoeff),it);
+ te_err_eff = gdat_data.te.error_bar(idatate(irhoeff),it);
+ % they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
+ ij=find(teeff./te_err_eff>10.);
+ if ~isempty(ij); te_err_eff(ij) = teeff(ij)./0.1; end
+ %
+ teeff = [teeff(1); teeff];
+ te_err_eff = [1e4; te_err_eff];
+ if gdat_data.gdat_params.fit_te_edge < 0
+ [gdat_data.fit.te.data(:,it), gdat_data.fit.te.d_over_drhopolnorm(:,it)] = ...
+ interpos(rhoeff,teeff,rhopolfit,fit_tension.te,[1 0],[0 0],te_err_eff);
+ else
+ [gdat_data.fit.te.data(:,it), gdat_data.fit.te.d_over_drhopolnorm(:,it)] = ...
+ interpos(rhoeff,teeff,rhopolfit,fit_tension.te,[1 2],[0 gdat_data.gdat_params.fit_te_edge],te_err_eff);
+ end
+ % impose minimum value
+ ij=find(gdat_data.fit.te.data(:,it)<fit_te_min);
+ if ~isempty(ij); gdat_data.fit.te.data(ij,it) = fit_te_min; end
end
- end
- end
- if any(strfind(data_request_eff,'ne'))
- idatane = find(gdat_data.ne.data(:,it)>1 & gdat_data.grids_1d.rhopolnorm(:,it)<=1.05);
- if length(idatane)>0
- gdat_data.fit.raw.ne.rhopolnorm(idatane,it) = gdat_data.grids_1d.rhopolnorm(idatane,it);
- gdat_data.fit.raw.ne.data(idatane,it) = gdat_data.ne.data(idatane,it);
- gdat_data.fit.raw.ne.error_bar(idatane,it) = gdat_data.ne.error_bar(idatane,it);
- [rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhopolnorm(idatane,it));
- rhoeff = [0; rhoeff];
- % they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
- neeff = gdat_data.ne.data(idatane(irhoeff),it);
- ne_err_eff = gdat_data.ne.error_bar(idatane(irhoeff),it);
- ij=find(neeff./ne_err_eff>10.);
- if ~isempty(ij); ne_err_eff(ij) = neeff(ij)./0.1; end
- %
- neeff = [neeff(1); neeff];
- ne_err_eff = [1e21; ne_err_eff];
- if gdat_data.gdat_params.fit_ne_edge < 0
- [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhopolnorm(:,it)] = interpos(rhoeff,neeff,rhopolfit,fit_tension.ne,[1 0],[0 ...
+ elseif strcmp(data_request_effi{i},'ne')
+ idatane = find(gdat_data.ne.data(indices_ok,it)>1 & gdat_data.grids_1d.rhopolnorm(indices_ok,it)<=1.05);
+ idatane = indices_ok(idatane);
+ if length(idatane)>0
+ gdat_data.fit.raw.ne.rhopolnorm(idatane,it) = gdat_data.grids_1d.rhopolnorm(idatane,it);
+ gdat_data.fit.raw.ne.data(idatane,it) = gdat_data.ne.data(idatane,it);
+ gdat_data.fit.raw.ne.error_bar(idatane,it) = gdat_data.ne.error_bar(idatane,it);
+ [rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhopolnorm(idatane,it));
+ rhoeff = [0; rhoeff];
+ % they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
+ neeff = gdat_data.ne.data(idatane(irhoeff),it);
+ ne_err_eff = gdat_data.ne.error_bar(idatane(irhoeff),it);
+ ij=find(neeff./ne_err_eff>10.);
+ if ~isempty(ij); ne_err_eff(ij) = neeff(ij)./0.1; end
+ %
+ neeff = [neeff(1); neeff];
+ ne_err_eff = [1e21; ne_err_eff];
+ if gdat_data.gdat_params.fit_ne_edge < 0
+ [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.d_over_drhopolnorm(:,it)] = interpos(rhoeff,neeff,rhopolfit,fit_tension.ne,[1 0],[0 ...
0],ne_err_eff);
- else
- ij_in_1 = find(rhoeff<1);
- [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhopolnorm(:,it)] = ...
- interpos([rhoeff(ij_in_1); 1.],[neeff(ij_in_1); gdat_data.gdat_params.fit_ne_edge], ...
+ else
+ ij_in_1 = find(rhoeff<1);
+ [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.d_over_drhopolnorm(:,it)] = ...
+ interpos([rhoeff(ij_in_1); 1.],[neeff(ij_in_1); gdat_data.gdat_params.fit_ne_edge], ...
rhopolfit,fit_tension.ne,[1 2],[0 gdat_data.gdat_params.fit_ne_edge],[ne_err_eff(ij_in_1);ne_err_eff(1)]);
+ end
+ % impose minimum value
+ ij=find(gdat_data.fit.ne.data(:,it)<fit_ne_min);
+ if ~isempty(ij); gdat_data.fit.ne.data(ij,it) = fit_ne_min; end
end
- end
- end
+ end
+ end
end
end
diff --git a/crpptbx/JET/jet_help_parameters.m b/crpptbx/JET/jet_help_parameters.m
index 0635755a..9b787201 100644
--- a/crpptbx/JET/jet_help_parameters.m
+++ b/crpptbx/JET/jet_help_parameters.m
@@ -38,6 +38,7 @@ help_struct_all.cocos = ['cocos value desired in output, uses eqdsk_cocos_transf
help_struct_all.fit = '0, no fit profiles, 1 (default) if fit profiles desired as well, relevant for _rho profiles. See also fit_type';
help_struct_all.fit_ne_edge = '-1, no edge value imposed, >0, value at boundary for all times';
help_struct_all.fit_te_edge = '-1, no edge value imposed, >0, value at boundary for all times';
+help_struct_all.fit_mask = '-1, no mask, [i1 i2 ...] index of channels to mask, plot(gdat_data.data) to check indices)';
% $$$ help_struct_all.fit_type = 'provenance of fitted profiles ''conf'' (default) from conf nodes, ''avg'' or ''local'' for resp. proffit: nodes';
help_struct_all.source = 'sxr: ''H'' (default), camera name ''T'', ''V'' ; for Vloop: ''ppf''(default),''jpf''';
help_struct_all.camera = ['[] (default, all), [i1 i2 ...] chord nbs ([1 3 5] if only chords 1, 3 and 5 are desired), ''central'' uses H10'];
diff --git a/crpptbx/JET/jet_requests_mapping.m b/crpptbx/JET/jet_requests_mapping.m
index f76d7674..d78e5600 100644
--- a/crpptbx/JET/jet_requests_mapping.m
+++ b/crpptbx/JET/jet_requests_mapping.m
@@ -188,6 +188,10 @@ switch lower(data_request)
mapping.timedim = 2;
mapping.label = 'ne';
mapping.method = 'switchcase';
+ case 'nete'
+ mapping.timedim = 2;
+ mapping.label = 'ne and Te'; % from chain2 or hrts with interpos fits
+ mapping.method = 'switchcase';
case 'nete_rho'
mapping.timedim = 2;
mapping.label = 'ne and Te';
--
GitLab