diff --git a/matlab/D3D/d3d_requests_mapping.m b/matlab/D3D/d3d_requests_mapping.m
index f3789d1c9ea3c580ad858d56dfa15119a416dda2..84879d6ffc32a7347a6dabd35dd18df55c10334a 100644
--- a/matlab/D3D/d3d_requests_mapping.m
+++ b/matlab/D3D/d3d_requests_mapping.m
@@ -109,6 +109,11 @@ switch lower(data_request)
mapping.gdat_timedim = 2;
mapping.method = 'switchcase';
mapping.expression = '';
+ case {'gas','gasa'}
+ mapping.timedim = 1;
+ mapping.label = 'gasa';
+ mapping.method = 'signal';
+ mapping.expression = [{'D3D'},{'\gasa'}];
case 'halpha'
mapping.timedim = 1;
mapping.label = 'Halpha';
@@ -194,19 +199,26 @@ switch lower(data_request)
mapping.timedim = 1;
mapping.label = 'line integrated el. density';
mapping.method = 'signal';
- mapping.expression = {'d3d','\d3d::top.electrons.bci:denv1'};
+ mapping.expression = {'d3d','\d3d::top.electrons.bci:denr0'};
case {'neints'}
mapping.timedim = 1;
- mapping.label = {'denv1','denv2','denv3'};
+ mapping.label = {'denr0','denv1','denv2','denv3'};
mapping.method = 'expression';
- mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''\d3d::top.electrons.bci:denv1'';' ...
+ mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''\d3d::top.electrons.bci:denr0'';' ...
'gdat_tmp=gdat_d3d(shot,params_eff);gdat_tmp.denv1=gdat_tmp;' ...
+ 'params_eff.data_request=''\d3d::top.electrons.bci:denv1'';gdat_denv1=gdat_d3d(shot,params_eff);' ...
+ 'gdat_tmp.denv1=gdat_denv1;gdat_tmp.data(:,2)=gdat_denv1.data;' ...
'params_eff.data_request=''\d3d::top.electrons.bci:denv2'';gdat_denv2=gdat_d3d(shot,params_eff);' ...
- 'gdat_tmp.denv2=gdat_denv2;gdat_tmp.data(:,2)=gdat_denv2.data;' ...
+ 'gdat_tmp.denv2=gdat_denv2;gdat_tmp.data(:,3)=gdat_denv2.data;' ...
'params_eff.data_request=''\d3d::top.electrons.bci:denv3'';gdat_denv3=gdat_d3d(shot,params_eff);' ...
- 'gdat_tmp.denv3=gdat_denv3;gdat_tmp.data(:,3)=gdat_denv3.data;' ...
- 'gdat_tmp.dim{2} = [1:3];gdat_tmp.dimunits{2} = {''denv1'',''denv2'',''denv3''};' ...
+ 'gdat_tmp.denv3=gdat_denv3;gdat_tmp.data(:,4)=gdat_denv3.data;' ...
+ 'gdat_tmp.dim{2} = [1:4];gdat_tmp.dimunits{2} = {''denr0'',''denv1'',''denv2'',''denv3''};' ...
';gdat_tmp.dimunits{1}=''s'';'];
+ case {'neint_ts'}
+ % assumes ne_rho in source input
+ mapping.timedim = 1;
+ mapping.label = 'neint from TS';
+ mapping.method = 'switchcase';
case 'nel'
mapping.timedim = 1;
mapping.label = 'NEBAR\_R0';
diff --git a/matlab/D3D/gdat_d3d.m b/matlab/D3D/gdat_d3d.m
index 89596d87d6e5bbdd0dfa8bc6281957440e43fe17..c3b5c90df0453e60fb371d841279af6b0afd3fe1 100644
--- a/matlab/D3D/gdat_d3d.m
+++ b/matlab/D3D/gdat_d3d.m
@@ -738,7 +738,13 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
gdat_data.psi_axis(it)= equil_all_t.gdata(it).ssimag;
gdat_data.psi_lcfs(it)= equil_all_t.gdata(it).ssibry;
gdat_data.rhopolnorm(:,it) = sqrt(abs((gdat_data.psi(:,it)-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it))));
+ gdat_data.rhopolnorm(1,it) = 0.;
gdat_data.rhopolnorm(end,it) = 1.;
+ [~,~,~,int_q]=interpos(gdat_data.psi(:,it),gdat_data.qvalue(:,it));
+ gdat_data.rhotornorm(:,it) = sqrt(int_q./int_q(end));
+ gdat_data.rhotornorm(1,it) = 0.;
+ gdat_data.rhotornorm(end,it) = 1.;
+ gdat_data.phi_intq(it) = int_q(end);
gdat_data.vol(:,it)=gdat_data.rhovolnorm(:,it).^2 .* eq.results.aeqdsk.volume(it);
gdat_data.Rmesh(:,it) = eq.results.geqdsk.r;
gdat_data.Zmesh(:,it) = eq.results.geqdsk.z;
@@ -786,6 +792,20 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ne','te'}
+ % Tree: ELECTRONS
+ %
+ % TOP.TS.BLESSED.CORE
+ % TOP.TS.BLESSED.DIVERTOR
+ % TOP.TS.BLESSED.TANGENTIAL
+ % TEMP
+ % DENSITY
+ % TEMP_E
+ % DENSITY_E
+ % Z
+ % R
+ % DCDATA
+ % PLDATA
+ % PLERROR
exp_name_eff = gdat_data.gdat_params.exp_name;
% ne or Te from Thomson data on raw z mesh vs (z,t)
if strcmp(data_request_eff,'te')
@@ -793,70 +813,220 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
else
node_leaf = 'DENSITY';
end
- % core
- nodenameeff = ['TS.BLESSED.CORE:' node_leaf];
- a = gdat_d3d(shot,{'electrons',nodenameeff});
- if isempty(a.data) || isempty(a.t)
- if gdat_params.nverbose>=3;
- a
- disp(['with data_request= ' data_request_eff])
- end
+ %
+ extra_source = {'core', 'tangential', 'divertor'};
+ gdat_data.with_data = {};
+ a_ip = gdat_d3d(shot,'ip');
+ if isempty(a_ip.t)
+ warning('no Ip?');
return
+ else
+ t_range = [a_ip.t(1) a_ip.t(end)];
+ end
+ for i=1:numel(extra_source)
+ main_source = ['TS.BLESSED.' extra_source{i} ':'];
+ nodenameeff = [main_source node_leaf];
+ a = gdat_d3d(shot,{'electrons',nodenameeff});
+ if isempty(a.data) || isempty(a.t)
+ if gdat_params.nverbose>=3;
+ a
+ disp(['with data_request= ' data_request_eff])
+ end
+ gdat_data.(extra_source{i}) = a;
+ else
+ ij=[a.t<t_range(1) | a.t>t_range(2)];
+ a.data(ij,:) = NaN; % not yet transposed
+ gdat_data.with_data{end+1} = extra_source{i};
+ gdat_data.(extra_source{i}).data = a.data';
+ gdat_data.(extra_source{i}).t = a.t;
+ gdat_data.(extra_source{i}).x = a.x;
+ gdat_data.(extra_source{i}).gdat_params.tree = a.gdat_params.tree;
+ if any(strcmp(fieldnames(a),'units'))
+ gdat_data.(extra_source{i}).units=a.units;
+ end
+ gdat_data.(extra_source{i}).r = mdsdata([main_source 'r']);
+ gdat_data.(extra_source{i}).z = mdsdata([main_source 'z']);
+ gdat_data.(extra_source{i}).error_bar = mdsdata([nodenameeff '_e'])';
+ gdat_data.(extra_source{i}).data_fullpath=[data_request_eff 'from electrons ' nodenameeff];
+ gdat_data.(extra_source{i}).dim=[{gdat_data.x};{gdat_data.t}];
+ if strcmp(extra_source{i},'tangential')
+ gdat_data.(extra_source{i}).dimunits=[{'R [m]'}; {'time [s]'}];
+ else
+ gdat_data.(extra_source{i}).dimunits=[{'Z [m]'}; {'time [s]'}];
+ end
+ end
end
- gdat_data.t = a.t;
- gdat_data.x = a.x;
- gdat_data.gdat_params.tree = a.gdat_params.tree;
- gdat_data.gdat_params.equil = a.gdat_params.equil;
- if any(strcmp(fieldnames(a),'units'))
- gdat_data.units=a.units;
- end
- gdat_data.core.data = a.data';
- gdat_data.core.t = a.t;
- gdat_data.core.x = a.x;
-% $$$ gdat_data.(lower(node_child_nameeff)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff)).data,gdat_data.(lower(node_child_nameeff)).data-alow.value);
+ %
+ % put core in main
+ i = 1;
+ gdat_data.data = gdat_data.(extra_source{i}).data;
+ gdat_data.error_bar = gdat_data.(extra_source{i}).error_bar;
+ gdat_data.t = gdat_data.(extra_source{i}).t;
+ gdat_data.x = gdat_data.(extra_source{i}).x;
+ gdat_data.r = gdat_data.(extra_source{i}).r;
+ gdat_data.z = gdat_data.(extra_source{i}).z;
+ gdat_data.dim = gdat_data.(extra_source{i}).dim;
+ gdat_data.dimunits = gdat_data.(extra_source{i}).dimunits;
+ gdat_data.units = gdat_data.(extra_source{i}).units;
+ gdat_data.data_fullpath = gdat_data.(extra_source{i}).data_fullpath;
+ gdat_data.(data_request_eff) = gdat_data;
+ %
-% $$$ [a,error_status]=rdaD3D_eff(shot,'VTA','R_core',exp_name_eff);
-% $$$ gdat_data.(lower(node_child_nameeff)).r = repmat(a.data,size(gdat_data.(lower(node_child_nameeff)).data,1),1);
-% $$$ [a,error_status]=rdaD3D_eff(shot,'VTA','Z_core',exp_name_eff);
- gdat_data.z = gdat_data.x;
- gdat_data.data_fullpath=[data_request_eff 'from ' nodenameeff];
- nb_core = numel(gdat_data.x);
- gdat_data.data = a.data'; % to get time as 2nd dim
- % gdat_data.data(1:nb_core,:) = gdat_data.data(1:nb_core,:);
- gdat_data.dim=[{gdat_data.z};{gdat_data.t}];
- gdat_data.dimunits=[{'Z [m]'}; {'time [s]'}];
-% $$$ gdat_data.error_bar(1:nb_core,:) = gdat_data.error_bar(1:nb_core,:);
- return
- % add edge part
- nodenameeff_e = [upper(data_request_eff(1)) 'e_e'];
- node_child_nameeff_e = [upper(data_request_eff(1)) 'e_edge'];
- [a,error_status]=rdaD3D_eff(shot,'VTA',nodenameeff_e,exp_name_eff);
- gdat_data.(lower(node_child_nameeff_e)).data = a.data;
- ij_edge_notok = find(a.data>5e21);
- gdat_data.(lower(node_child_nameeff_e)).data(ij_edge_notok) = NaN;
- gdat_data.(lower(node_child_nameeff_e)).t = a.t;
- if ~isempty(a.data)
- [a,error_status]=rdaD3D_eff(shot,'VTA','R_edge',exp_name_eff);
- gdat_data.(lower(node_child_nameeff_e)).r = repmat(a.data,size(gdat_data.(lower(node_child_nameeff_e)).data,1),1);
- [a,error_status]=rdaD3D_eff(shot,'VTA','Z_edge',exp_name_eff);
- gdat_data.(lower(node_child_nameeff_e)).z = repmat(a.data',1,size(gdat_data.(lower(node_child_nameeff_e)).data,2));
- nb_edge = size(gdat_data.(lower(node_child_nameeff_e)).z,1);
- iaaa=iround_os(gdat_data.(lower(node_child_nameeff_e)).t,gdat_data.(lower(node_child_nameeff)).t);
- gdat_data.data(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).data(1:nb_edge,iaaa);
- gdat_data.dim{1}(nb_core+1:nb_core+nb_edge,:)=gdat_data.(lower(node_child_nameeff_e)).z(1:nb_edge,iaaa);
- [alow,e]=rdaD3D_eff(shot,'VTA',[upper(data_request_eff(1)) 'elow_e'],exp_name_eff);
- [aup,e]=rdaD3D_eff(shot,'VTA',[upper(data_request_eff(1)) 'eupp_e'],exp_name_eff);
- gdat_data.(lower(node_child_nameeff_e)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff_e)).data, ...
- gdat_data.(lower(node_child_nameeff_e)).data-alow.value);
- gdat_data.error_bar(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).error_bar(1:nb_edge,iaaa);
- else
- gdat_data.(lower(node_child_nameeff_e)).data = [];
- gdat_data.(lower(node_child_nameeff_e)).t = [];
- gdat_data.(lower(node_child_nameeff_e)).error_bar = [];
- gdat_data.(lower(node_child_nameeff_e)).r = [];
- gdat_data.(lower(node_child_nameeff_e)).z = [];
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'neint_ts'}
+ if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || ~isstruct(gdat_data.gdat_params.source)
+ error('requires ''ne_rho'' gdat structure as input in ''source'' option');
+ end
+ ne_rho = gdat_data.gdat_params.source;
+ % Computes R0 and V1 line integrated signal based on ne fit (so inside LCFS), provide also R(psi_ts) and Z(psi_ts) on these lines
+ z_at_r0line = 0.;
+ r_at_v1line = 1.48;
+ % finds xx 1d dependence for each line within d3d wall for each time in equil
+ nb_xx_points = 257;
+ for it_equil=1:numel(ne_rho.equil.t)
+ % R0 line
+ gdat_data.r0.y = z_at_r0line;
+ gdat_data.r0.x(:,it_equil) = linspace(ne_rho.equil.Rmesh(1,it_equil),ne_rho.equil.Rmesh(end,it_equil),nb_xx_points);
+ [psi_xx_out]=interpos2Dcartesian(ne_rho.equil.Rmesh(:,it_equil),ne_rho.equil.Zmesh(:,it_equil), ...
+ ne_rho.equil.psi2D(:,:,it_equil),gdat_data.r0.x(:,it_equil),z_at_r0line.*ones(size(gdat_data.r0.x(:,it_equil))));
+ gdat_data.r0.psi(:,it_equil) = psi_xx_out;
+ gdat_data.r0.rhopolnorm(:,it_equil) = sqrt((psi_xx_out-ne_rho.equil.psi_axis(it_equil)) / ...
+ (ne_rho.equil.psi_lcfs(it_equil)-ne_rho.equil.psi_axis(it_equil)));
+ % find line within lcfs
+ zdist = (psi_xx_out(1:end-1)-ne_rho.equil.psi_lcfs(it_equil)).*(psi_xx_out(2:end)-ne_rho.equil.psi_lcfs(it_equil));
+ ij=find(zdist<=0);
+ % assume 2 crossing not near ends
+ for ii=1:2
+ r_cross_lcfs(ii) = interp1(psi_xx_out(ij(ii):ij(ii)+1),gdat_data.r0.x(ij(ii):ij(ii)+1,it_equil), ...
+ ne_rho.equil.psi_lcfs(it_equil));
+ end
+ gdat_data.r0.x_cross_lcfs(1:2,it_equil) = r_cross_lcfs;
+ gdat_data.r0.x_len_plasma(it_equil) = abs(diff(gdat_data.r0.x_cross_lcfs(:,it_equil)));
+ % find line within walls
+ zdist=(ne_rho.equil.zwall(1:end-1,it_equil)-z_at_r0line).*(ne_rho.equil.zwall(2:end,it_equil)-z_at_r0line);
+ ij=find(zdist<=0);
+ % assume 2 crossing
+ for ii=1:min(2,numel(ij))
+ if ij(ii)~=size(ne_rho.equil.zwall,1)
+ rcrosswall(ii) = interp1(ne_rho.equil.zwall(ij(ii):ij(ii)+1,it_equil),ne_rho.equil.rwall(ij(ii):ij(ii)+1,it_equil), ...
+ z_at_r0line);
+ else
+ rcrosswall(ii) = interp1(ne_rho.equil.zwall(ij(ii)-1:ij(ii),it_equil),ne_rho.equil.rwall(ij(ii)-1:ij(ii),it_equil), ...
+ z_at_r0line);
+ end
+ end
+ if numel(ij)==1
+ rcrosswall(2) = NaN;
+ end
+ gdat_data.r0.x_crosswall(1:2,it_equil) = rcrosswall;
+ gdat_data.r0.x_len_crosswall(it_equil) = abs(diff(gdat_data.r0.x_crosswall(:,it_equil)));
+ % V1 line
+ gdat_data.v1.y = r_at_v1line;
+ gdat_data.v1.x(:,it_equil) = linspace(ne_rho.equil.Zmesh(1,it_equil),ne_rho.equil.Zmesh(end,it_equil),nb_xx_points);
+ [psi_xx_out]=interpos2Dcartesian(ne_rho.equil.Rmesh(:,it_equil),ne_rho.equil.Zmesh(:,it_equil), ...
+ ne_rho.equil.psi2D(:,:,it_equil),r_at_v1line.*ones(size(gdat_data.v1.x(:,it_equil))),gdat_data.v1.x(:,it_equil));
+ gdat_data.v1.psi(:,it_equil) = psi_xx_out;
+ gdat_data.v1.rhopolnorm(:,it_equil) = sqrt((psi_xx_out-ne_rho.equil.psi_axis(it_equil)) / ...
+ (ne_rho.equil.psi_lcfs(it_equil)-ne_rho.equil.psi_axis(it_equil)));
+ % find line within lcfs
+ zdist = (psi_xx_out(1:end-1)-ne_rho.equil.psi_lcfs(it_equil)).*(psi_xx_out(2:end)-ne_rho.equil.psi_lcfs(it_equil));
+ ij=find(zdist<=0);
+ % may have more than 2 crossing if private flux
+ if numel(ij)>2
+ zz_ij=abs(gdat_data.v1.x(ij,it_equil));
+ [dok,ijok] = sort(zz_ij);
+ else
+ ijok = [1:2];
+ end
+ for ii=1:2
+ r_cross_lcfs(ii) = interp1(psi_xx_out(ij(ijok(ii)):ij(ijok(ii))+1),gdat_data.v1.x(ij(ijok(ii)):ij(ijok(ii))+1,it_equil), ...
+ ne_rho.equil.psi_lcfs(it_equil));
+ end
+ gdat_data.v1.x_cross_lcfs(1:2,it_equil) = r_cross_lcfs;
+ gdat_data.v1.x_len_plasma(it_equil) = abs(diff(gdat_data.v1.x_cross_lcfs(:,it_equil)));
+ % walls
+ zdist=(ne_rho.equil.rwall(1:end-1,it_equil)-r_at_v1line).*(ne_rho.equil.rwall(2:end,it_equil)-r_at_v1line);
+ ij=find(zdist<=0 & ne_rho.equil.zwall(1:end-1,it_equil)<=0);
+ if numel(ij)>1
+ rho_ij=sqrt((ne_rho.equil.rwall(ij,it_equil)-r_at_v1line).^2 + (ne_rho.equil.zwall(ij,it_equil)-0.).^2);
+ [dok,ijok] = sort(rho_ij);
+ zcrosswall(1) = interp1(ne_rho.equil.rwall(ij(ijok):ij(ijok)+1,it_equil),ne_rho.equil.zwall(ij(ijok):ij(ijok)+1,it_equil), ...
+ r_at_v1line);
+ elseif numel(ij)==1
+ zcrosswall(1) = interp1(ne_rho.equil.rwall(ij:ij+1,it_equil),ne_rho.equil.zwall(ij:ij+1,it_equil), ...
+ r_at_v1line);
+ else
+ zcrosswall(1) = NaN;
+ end
+ ij=find(zdist<=0 & ne_rho.equil.zwall(1:end-1,it_equil)>=0);
+ if numel(ij)>1
+ rho_ij=sqrt((ne_rho.equil.rwall(ij,it_equil)-r_at_v1line).^2 + (ne_rho.equil.zwall(ij,it_equil)-0.).^2);
+ [dok,ijok] = sort(rho_ij);
+ zcrosswall(2) = interp1(ne_rho.equil.rwall(ij(ijok):ij(ijok)+1,it_equil),ne_rho.equil.zwall(ij(ijok):ij(ijok)+1,it_equil), ...
+ r_at_v1line);
+ elseif numel(ij)==1
+ zcrosswall(2) = interp1(ne_rho.equil.rwall(ij:ij+1,it_equil),ne_rho.equil.zwall(ij:ij+1,it_equil), ...
+ r_at_v1line);
+ else
+ zcrosswall(2) = NaN;
+ end
+ gdat_data.v1.x_crosswall(1:2,it_equil) = zcrosswall;
+ gdat_data.v1.x_len_crosswall(it_equil) = abs(diff(gdat_data.v1.x_crosswall(:,it_equil)));
+ end
+ for isub={'r0','v1'}
+ gdat_data.(isub{1}).t = ne_rho.fit.t;
+ gdat_data.(isub{1}).dimunits{1} = 'time [s]';
+ gdat_data.(isub{1}).dim{1} = gdat_data.(isub{1}).t;
+ gdat_data.(isub{1}).units = 'm^{-2}';
+ gdat_data.(isub{1}).label = ['ne int one line ' upper(isub{1})];
+ gdat_data.(isub{1}).data_fullpath = ['from ne_rho fit to Thomson line-integral along line ' upper(isub{1})];
+ end
+ for it=1:numel(ne_rho.fit.t)
+ it_equil = iround_os(ne_rho.equil.t,ne_rho.fit.t(it));
+ % r0
+ gdat_data.r0.it_equil(it) = it_equil;
+ gdat_data.r0.nefit(:,it) = interpos(-63,ne_rho.fit.rhopolnorm(:,it),ne_rho.fit.ne.data(:,it),gdat_data.r0.rhopolnorm(:,it_equil), ...
+ -0.1,[1 0],[0 0]);
+ ij = find((gdat_data.r0.x(:,it_equil)-gdat_data.r0.x_crosswall(1,it_equil)).*(gdat_data.r0.x(:,it_equil)-gdat_data.r0.x_crosswall(2,it_equil)) > 0);
+ gdat_data.r0.nefit(ij,it) = NaN;
+ ij=find(~isnan(gdat_data.r0.nefit(:,it)));
+ gdat_data.r0.data(it) = trapz(gdat_data.r0.x(ij,it_equil),gdat_data.r0.nefit(ij,it));
+ non_vacuum = gdat_data.r0.nefit(:,it);
+ non_vacuum_min = min(non_vacuum(non_vacuum>0));
+ non_vacuum(non_vacuum==0) = non_vacuum_min;
+ gdat_data.r0.neint_ts_vacuum_min(it) = trapz(gdat_data.r0.x(ij,it_equil),non_vacuum(ij));
+ non_vacuum_110 = gdat_data.r0.nefit(:,it);
+ ij = find(gdat_data.r0.rhopolnorm(:,it_equil)>1 & gdat_data.r0.rhopolnorm(:,it_equil)<=1.10 & ~isnan(non_vacuum_110));
+ non_vacuum_110(ij) = non_vacuum_min;
+ ij=find(~isnan(non_vacuum_110));
+ gdat_data.r0.neint_ts_vacuum_rhopol110(it) = trapz(gdat_data.r0.x(ij,it_equil),non_vacuum_110(ij));
+ % v1
+ gdat_data.v1.it_equil(it) = it_equil;
+ gdat_data.v1.nefit(:,it) = interpos(-63,ne_rho.fit.rhopolnorm(:,it),ne_rho.fit.ne.data(:,it),gdat_data.v1.rhopolnorm(:,it_equil), ...
+ -0.1,[1 0],[0 0]);
+ ij = find((gdat_data.v1.x(:,it_equil)-gdat_data.v1.x_crosswall(1,it_equil)).*(gdat_data.v1.x(:,it_equil)-gdat_data.v1.x_crosswall(2,it_equil)) > 0);
+ gdat_data.v1.nefit(ij,it) = NaN;
+ ij=find(~isnan(gdat_data.v1.nefit(:,it)));
+ gdat_data.v1.data(it) = trapz(gdat_data.v1.x(ij,it_equil),gdat_data.v1.nefit(ij,it));
+ non_vacuum = gdat_data.v1.nefit(:,it);
+ non_vacuum_min = min(non_vacuum(non_vacuum>0));
+ non_vacuum(non_vacuum==0) = non_vacuum_min;
+ gdat_data.v1.neint_ts_vacuum_min(it) = trapz(gdat_data.v1.x(ij,it_equil),non_vacuum(ij));
+ non_vacuum_110 = gdat_data.v1.nefit(:,it);
+ ij = find(gdat_data.v1.rhopolnorm(:,it_equil)>1 & gdat_data.v1.rhopolnorm(:,it_equil)<=1.1 & ~isnan(non_vacuum_110));
+ non_vacuum_110(ij) = non_vacuum_min;
+ ij=find(~isnan(non_vacuum_110));
+ gdat_data.v1.neint_ts_vacuum_rhopol110(it) = trapz(gdat_data.v1.x(ij,it_equil),non_vacuum_110(ij));
+ end
+ gdat_data.r0.nel = gdat_data.r0.data./gdat_data.r0.x_len_plasma(gdat_data.r0.it_equil);
+ gdat_data.v1.nel = gdat_data.v1.data./gdat_data.v1.x_len_plasma(gdat_data.v1.it_equil);
+
+ fields_to_copy = {'data','units','dim','dimunits','t','x','data_fullpath','label','help'};
+ for i=1:length(fields_to_copy)
+ if isfield(gdat_data.r0,fields_to_copy{i})
+ gdat_data.(fields_to_copy{i}) = gdat_data.r0.(fields_to_copy{i});
+ end
end
- gdat_data.x=gdat_data.dim{1};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ne_rho', 'te_rho', 'nete_rho'}
@@ -869,6 +1039,11 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
[gdat_data,params_kin,error_status]=gdat_d3d(shot,params_eff);
gdat_data.gdat_params.data_request=data_request_eff;
gdat_data.gdat_request=data_request_eff;
+ if strcmp(lower(data_request_eff(1:4)),'nete')
+ params_eff.data_request=data_request_eff(3:4); % start with ne if nete_rho
+ [aate,params_kin,error_status]=gdat_d3d(shot,params_eff);
+ gdat_data.te = aate.te;
+ end
if error_status>0
if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
return
@@ -881,117 +1056,61 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
return
end
gdat_data.gdat_params.equil = params_equil.equil;
- %gdat_data.equil = equil;
+ gdat_data.equil = equil;
- % core
- node_child_nameeff = [upper(data_request_eff(1)) 'e_core'];
- inb_chord_core=size(gdat_data.(lower(node_child_nameeff)).r,1);
- inb_time_core=size(gdat_data.(lower(node_child_nameeff)).r,2);
- psi_out_core = NaN*ones(inb_chord_core,inb_time_core);
- rhopolnorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
- rhotornorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
- rhovolnorm_out_core = NaN*ones(inb_chord_core,inb_time_core);
- % edge
- node_child_nameeff_e = [upper(data_request_eff(1)) 'e_edge'];
- inb_chord_edge=size(gdat_data.(lower(node_child_nameeff_e)).r,1);
- inb_time_edge=size(gdat_data.(lower(node_child_nameeff_e)).r,2);
- psi_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
- rhopolnorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
- rhotornorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
- rhovolnorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
- % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
- time_equil=[min(gdat_data.(lower(node_child_nameeff)).t(1)-0.1,equil.t(1)-0.1) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) max(equil.t(end)+0.1,gdat_data.(lower(node_child_nameeff)).t(end)+0.1)];
- for itequil=1:length(time_equil)-1
- rr=equil.Rmesh(:,itequil);
- zz=equil.Zmesh(:,itequil);
- psirz_in = equil.psi2D(:,:,itequil);
- it_core_inequil = find(gdat_data.(lower(node_child_nameeff)).t>time_equil(itequil) & gdat_data.(lower(node_child_nameeff)).t<=time_equil(itequil+1));
- if ~isempty(it_core_inequil)
- rout_core=gdat_data.(lower(node_child_nameeff)).r(:,it_core_inequil);
- zout_core=gdat_data.(lower(node_child_nameeff)).z(:,it_core_inequil);
- psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core);
- psi_out_core(:,it_core_inequil) = reshape(psi_at_routzout,inb_chord_core,length(it_core_inequil));
- rhopolnorm_out_core(:,it_core_inequil) = sqrt((psi_out_core(:,it_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
- for it_cx=1:length(it_core_inequil)
- rhotornorm_out_core(:,it_core_inequil(it_cx)) = ...
- interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
- rhovolnorm_out_core(:,it_core_inequil(it_cx)) = ...
- interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
- end
- end
- % edge
- it_edge_inequil = find(gdat_data.(lower(node_child_nameeff_e)).t>time_equil(itequil) & gdat_data.(lower(node_child_nameeff_e)).t<=time_equil(itequil+1));
- if ~isempty(it_edge_inequil)
- rout_edge=gdat_data.(lower(node_child_nameeff_e)).r(:,it_edge_inequil);
- zout_edge=gdat_data.(lower(node_child_nameeff_e)).z(:,it_edge_inequil);
- psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_edge,zout_edge);
- psi_out_edge(:,it_edge_inequil) = reshape(psi_at_routzout,inb_chord_edge,length(it_edge_inequil));
- rhopolnorm_out_edge(:,it_edge_inequil) = sqrt((psi_out_edge(:,it_edge_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
- for it_cx=1:length(it_edge_inequil)
- rhotornorm_out_edge(:,it_edge_inequil(it_cx)) = ...
- interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
- rhovolnorm_out_edge(:,it_edge_inequil(it_cx)) = ...
- interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
+ % sources
+ % assume, time and grid same for ne and Te for each extra_source
+ extra_source = {'core', 'tangential', 'divertor'};
+ for i=1:numel(extra_source)
+ node_child_nameeff = extra_source{i};
+ if ~isempty(gdat_data.(extra_source{i}).data) && ~ischar(gdat_data.(extra_source{i}).data)
+ inb_chord_source=size(gdat_data.(extra_source{i}).data,1);
+ inb_time_source=size(gdat_data.(extra_source{i}).data,2);
+ psi_out_source = NaN*ones(inb_chord_source,inb_time_source);
+ rhopolnorm_out_source = NaN*ones(inb_chord_source,inb_time_source);
+ rhotornorm_out_source = NaN*ones(inb_chord_source,inb_time_source);
+ rhovolnorm_out_source = NaN*ones(inb_chord_source,inb_time_source);
+ % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
+ time_equil=[min(gdat_data.(extra_source{i}).t(1)-0.1,equil.t(1)-0.1); 0.5.*(equil.t(1:end-1)+equil.t(2:end)) ;max(equil.t(end)+0.1,gdat_data.(extra_source{i}).t(end)+0.1)];
+ clear psi_out_core rhopolnorm_out_core rhotornorm_out_core rhovolnorm_out_core
+ for itequil=1:numel(time_equil)-1
+ rr=equil.Rmesh(:,itequil);
+ zz=equil.Zmesh(:,itequil);
+ psirz_in = equil.psi2D(:,:,itequil);
+ it_core_inequil = find(gdat_data.(extra_source{i}).t>time_equil(itequil) & gdat_data.(extra_source{i}).t<=time_equil(itequil+1));
+ if ~isempty(it_core_inequil)
+ rout_core=gdat_data.(extra_source{i}).r;
+ zout_core=gdat_data.(extra_source{i}).z;
+ psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core);
+ psi_out_core(:,it_core_inequil) = repmat(psi_at_routzout,1,numel(it_core_inequil));
+ rhopolnorm_out_core(:,it_core_inequil) = sqrt((psi_out_core(:,it_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
+ for it_cx=1:length(it_core_inequil)
+ rhotornorm_out_core(:,it_core_inequil(it_cx)) = ...
+ interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
+ rhovolnorm_out_core(:,it_core_inequil(it_cx)) = ...
+ interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
+ end
+ end
end
end
- end
- gdat_data.(lower(node_child_nameeff)).psi = psi_out_core;
- gdat_data.(lower(node_child_nameeff)).rhopolnorm = rhopolnorm_out_core;
- gdat_data.(lower(node_child_nameeff)).rhotornorm = rhotornorm_out_core;
- gdat_data.(lower(node_child_nameeff)).rhovolnorm = rhovolnorm_out_core;
- gdat_data.(lower(node_child_nameeff_e)).psi = psi_out_edge;
- gdat_data.(lower(node_child_nameeff_e)).rhopolnorm = rhopolnorm_out_edge;
- gdat_data.(lower(node_child_nameeff_e)).rhotornorm = rhotornorm_out_edge;
- gdat_data.(lower(node_child_nameeff_e)).rhovolnorm = rhovolnorm_out_edge;
- % put values of rhopolnorm for dim{1} by default, all radial mesh for combined core, edge in grids_1d
- gdat_data.x = gdat_data.(lower(node_child_nameeff)).rhopolnorm;
- iaaa=iround_os(gdat_data.(lower(node_child_nameeff_e)).t,gdat_data.(lower(node_child_nameeff)).t);
- gdat_data.x(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhopolnorm(:,iaaa);
- gdat_data.dim{1} = gdat_data.x;
- gdat_data.dimunits{1} = 'rhopolnorm';
- gdat_data.grids_1d.rhopolnorm = gdat_data.x;
- gdat_data.grids_1d.psi = gdat_data.(lower(node_child_nameeff)).psi;
- gdat_data.grids_1d.psi(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).psi(:,iaaa);
- gdat_data.grids_1d.rhotornorm = gdat_data.(lower(node_child_nameeff)).rhotornorm;
- gdat_data.grids_1d.rhotornorm(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhotornorm(:,iaaa);
- gdat_data.grids_1d.rhovolnorm = gdat_data.(lower(node_child_nameeff)).rhovolnorm;
- gdat_data.grids_1d.rhovolnorm(inb_chord_core+1:inb_chord_core+inb_chord_edge,:) = gdat_data.(lower(node_child_nameeff_e)).rhovolnorm(:,iaaa);
-
- gdat_data.data_fullpath = [gdat_data.data_fullpath ' projected on equilibrium ' gdat_data.gdat_params.equil];
+ gdat_data.(extra_source{i}).psi = psi_out_core;
+ gdat_data.(extra_source{i}).rhopolnorm = rhopolnorm_out_core;
+ gdat_data.(extra_source{i}).rhotornorm = rhotornorm_out_core;
+ gdat_data.(extra_source{i}).rhovolnorm = rhovolnorm_out_core;
+ % put values of rhopolnorm for dim{1} by default, all radial mesh for combined core, edge in grids_1d
+ if i==1
+ gdat_data.x = gdat_data.(extra_source{i}).rhopolnorm;
+ gdat_data.dim{1} = gdat_data.x;
+ gdat_data.dimunits{1} = 'rhopolnorm';
+ gdat_data.grids_1d.rhopolnorm = gdat_data.x;
+ gdat_data.grids_1d.psi = gdat_data.(extra_source{i}).psi;
+ gdat_data.grids_1d.rhotornorm = gdat_data.(extra_source{i}).rhotornorm;
+ gdat_data.grids_1d.rhovolnorm = gdat_data.(extra_source{i}).rhovolnorm;
- % if nete_rho, copy data as .ne, get .te and put pe=e ne Te in data:
- gdat_data.(data_request_eff(1:2)).data = gdat_data.data;
- gdat_data.(data_request_eff(1:2)).error_bar = gdat_data.error_bar;
- gdat_data.(data_request_eff(1:2)).units = gdat_data.units;
- gdat_data.(data_request_eff(1:2)).core = gdat_data.([data_request_eff(1:2) '_core']);
- gdat_data.(data_request_eff(1:2)).edge = gdat_data.([data_request_eff(1:2) '_edge']);
- gdat_data = rmfield(gdat_data,{[data_request_eff(1:2) '_core'],[data_request_eff(1:2) '_edge']});
- if strcmp(data_request_eff(1:4),'nete')
- params_eff.data_request=data_request_eff(3:4);
- [gdat_data_te,params_kin,error_status]=gdat_d3d(shot,params_eff);
- gdat_data.te.data = gdat_data_te.data;
- gdat_data.te.error_bar = gdat_data_te.error_bar;
- gdat_data.te.units = gdat_data_te.units;
- gdat_data.te.core = gdat_data_te.te_core;
- gdat_data.te.edge = gdat_data_te.te_edge;
- gdat_data.te.error_bar = gdat_data_te.error_bar;
- gdat_data.te.core.psi = gdat_data.ne.core.psi;
- gdat_data.te.core.rhopolnorm = gdat_data.ne.core.rhopolnorm;
- gdat_data.te.core.rhotornorm = gdat_data.ne.core.rhotornorm;
- gdat_data.te.core.rhovolnorm = gdat_data.ne.core.rhovolnorm;
- gdat_data.te.edge.psi = gdat_data.ne.edge.psi;
- gdat_data.te.edge.rhopolnorm = gdat_data.ne.edge.rhopolnorm;
- gdat_data.te.edge.rhotornorm = gdat_data.ne.edge.rhotornorm;
- gdat_data.te.edge.rhovolnorm = gdat_data.ne.edge.rhovolnorm;
- % put pe in main gdat_data
- gdat_data.data = 1.6022e-19.*gdat_data.ne.data.*gdat_data.te.data;
- 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);
- gdat_data.units='N/m^2; 1.6022e-19 ne Te';
- gdat_data.data_fullpath=['pe=1.6e-19*ne*Te in data, .ne, .te ' gdat_data.data_fullpath(3:end)];
- gdat_data.label = 'pe';
+ gdat_data.data_fullpath = [gdat_data.data_fullpath ' projected on equilibrium ' gdat_data.gdat_params.equil];
+ end
end
-
+ %
% defaults for fits, so user always gets std structure
gdat_data.fit.rhotornorm = []; % same for both ne and te
gdat_data.fit.rhopolnorm = [];
@@ -1053,7 +1172,7 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
gdat_data.fit.rhopolnorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.grids_1d.rhopolnorm(irhook(isort),it); 1],rhotornormfit,-0.1,[2 2],[0 1]);
if any(strfind(data_request_eff,'te'))
idatate = find(gdat_data.te.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05);
- if length(idatate)>0
+ if numel(idatate) >= 1
gdat_data.fit.raw.te.rhotornorm(idatate,it) = gdat_data.grids_1d.rhotornorm(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);
@@ -1072,7 +1191,7 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
end
if any(strfind(data_request_eff,'ne'))
idatane = find(gdat_data.ne.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05);
- if length(idatane)>0
+ if numel(idatane)>= 1
gdat_data.fit.raw.ne.rhotornorm(idatane,it) = gdat_data.grids_1d.rhotornorm(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);
@@ -1426,9 +1545,9 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
gdat_data.qvalue(end,it) = 2. * q_edge;
end
% $$$ if qpsi.value(ijok(1),1)<0
-% $$$ gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
+% $$$ gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
% $$$ else
-% $$$ gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
+% $$$ gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
% $$$ end
end
gdat_data.x = gdat_data.rhopolnorm;