Skip to content
Snippets Groups Projects
Commit 65b6b29b authored by Olivier Sauter's avatar Olivier Sauter
Browse files

add protection if part of the data is missing, so can do gdat('te_rho' with...

add protection if part of the data is missing, so can do gdat('te_rho' with only iter=1, rt-like, liuqe results without psitbx calculations for tcv

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@8484 d63d8f72-b253-0410-a779-e742ad2e26cf
parent ae19f86d
No related branches found
No related tags found
No related merge requests found
...@@ -422,7 +422,7 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi') ...@@ -422,7 +422,7 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi')
do_liuqem2liuqef = 1; do_liuqem2liuqef = 1;
end end
end end
eval_expr = ['tdi(''' mapping_for_tcv.expression ''');'] eval_expr = ['tdi(''' mapping_for_tcv.expression ''');'];
end end
% special cases for matching liuqe fortran and matlab % special cases for matching liuqe fortran and matlab
if liuqe_matlab == 0 && do_liuqem2liuqef==1 if liuqe_matlab == 0 && do_liuqem2liuqef==1
...@@ -1532,7 +1532,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ...@@ -1532,7 +1532,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
phi_tor.data = []; phi_tor.data = [];
phi_tor.units = 'Wb'; phi_tor.units = 'Wb';
end end
if ~any(any(isfinite(phi_tor.data))) if ~any(any(isfinite(phi_tor.data))) || ischar(phi_tor.data)
% no phi_tor, compute it from q profile % no phi_tor, compute it from q profile
q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']); q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']);
if liuqe_matlab==0 if liuqe_matlab==0
...@@ -1542,15 +1542,19 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ...@@ -1542,15 +1542,19 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
phi_tor.dimunits = q_rho.dimunits; phi_tor.dimunits = q_rho.dimunits;
params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE
psi_axis=gdat_tcv(shot,params_eff); psi_axis=gdat_tcv(shot,params_eff);
for it=1:size(q_rho.data,2) if isnumeric(q_rho.data)
ij=find(isfinite(q_rho.data(:,it))); for it=1:size(q_rho.data,2)
if ~isempty(ij) && numel(ij)>5 ij=find(isfinite(q_rho.data(:,it)));
[qfit,~,~,phi]=interpos(-43,q_rho.dim{1}(ij).^2,q_rho.data(ij,it),phi_tor.dim{1}.^2); if ~isempty(ij) && numel(ij)>5
% Phi=sigma_bp*sigma_rhothetaphi*(2pi)^(1-eBp)*int(q dpsi), dpsi=dpsi_norm * (psiedge-psaxis) [qfit,~,~,phi]=interpos(-43,q_rho.dim{1}(ij).^2,q_rho.data(ij,it),phi_tor.dim{1}.^2);
% cocos=17: sigma_Bp=-1,sigma_rhothetaphi=1, eBp=1, (psiedge-psaxis)=-psiaxis % Phi=sigma_bp*sigma_rhothetaphi*(2pi)^(1-eBp)*int(q dpsi), dpsi=dpsi_norm * (psiedge-psaxis)
phi = -1.* phi .* (-psi_axis.data(it)); % cocos=17: sigma_Bp=-1,sigma_rhothetaphi=1, eBp=1, (psiedge-psaxis)=-psiaxis
phi_tor.data(:,it) = phi'; phi = -1.* phi .* (-psi_axis.data(it));
phi_tor.data(:,it) = phi';
end
end end
else
phi_tor.data = [];
end end
elseif any(any(~isfinite(phi_tor.data))) elseif any(any(~isfinite(phi_tor.data)))
% there are some non-finite values, usually edge values, so replace them with integrated values % there are some non-finite values, usually edge values, so replace them with integrated values
...@@ -1588,8 +1592,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ...@@ -1588,8 +1592,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
gdat_data.dim = phi_tor.dim; gdat_data.dim = phi_tor.dim;
gdat_data.dimunits = phi_tor.dimunits; gdat_data.dimunits = phi_tor.dimunits;
gdat_data.units = phi_tor.units; gdat_data.units = phi_tor.units;
gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim}; if (length(gdat_data.dim)>=mapping_for_tcv.gdat_timedim); gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim}; end
gdat_data.x = gdat_data.dim{1}; if (length(gdat_data.dim)>0); gdat_data.x = gdat_data.dim{1}; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'psi_edge'} case {'psi_edge'}
...@@ -1711,32 +1715,37 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ...@@ -1711,32 +1715,37 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
b0=gdat_tcv([],params_eff); b0=gdat_tcv([],params_eff);
gdat_data.t = phi_tor.t; gdat_data.t = phi_tor.t;
ij=find(isfinite(b0.data)); ij=find(isfinite(b0.data));
gdat_data.b0 = interpos(b0.t(ij),b0.data(ij),gdat_data.t); if length(gdat_data.t) > 0
% always provide rhotor_edge so field always exists gdat_data.b0 = interpos(b0.t(ij),b0.data(ij),gdat_data.t);
gdat_data.rhotor_edge = sqrt(phi_tor.data(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)))'; % always provide rhotor_edge so field always exists
if strcmp(data_request_eff,'rhotor_edge') gdat_data.rhotor_edge = sqrt(phi_tor.data(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)))';
gdat_data.data = gdat_data.rhotor_edge; if strcmp(data_request_eff,'rhotor_edge')
gdat_data.dim = phi_tor.dim(2); % while there is a problem with \tcv_shot::top.results.equil... dim nodes gdat_data.data = gdat_data.rhotor_edge;
gdat_data.dimunits = phi_tor.dimunits{2}; gdat_data.dim = phi_tor.dim(2); % while there is a problem with \tcv_shot::top.results.equil... dim nodes
gdat_data.units = 'm'; gdat_data.dimunits = phi_tor.dimunits{2};
gdat_data.data_fullpath=['rhotor_edge=sqrt(Phi(edge)/pi/B0) from phi_tor']; gdat_data.units = 'm';
elseif strcmp(data_request_eff,'rhotor_norm') gdat_data.data_fullpath=['rhotor_edge=sqrt(Phi(edge)/pi/B0) from phi_tor'];
gdat_data.data = sqrt(phi_tor.data./(ones(size(phi_tor.data,1),1)*reshape(phi_tor.data(end,:),1,size(phi_tor.data,2)))); elseif strcmp(data_request_eff,'rhotor_norm')
gdat_data.dim = phi_tor.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes gdat_data.data = sqrt(phi_tor.data./(ones(size(phi_tor.data,1),1)*reshape(phi_tor.data(end,:),1,size(phi_tor.data,2))));
gdat_data.x=gdat_data.dim{1}; gdat_data.dim = phi_tor.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes
gdat_data.units = ''; gdat_data.x=gdat_data.dim{1};
gdat_data.dimunits = phi_tor.dimunits; gdat_data.units = '';
gdat_data.data_fullpath=['rhotor_norm=sqrt(Phi/Phi(edge)) from phi_tor']; gdat_data.dimunits = phi_tor.dimunits;
elseif strcmp(data_request_eff,'rhotor') gdat_data.data_fullpath=['rhotor_norm=sqrt(Phi/Phi(edge)) from phi_tor'];
gdat_data.data = sqrt(phi_tor.data./pi./(ones(size(phi_tor.data,1),1)*reshape(gdat_data.b0,1,numel(gdat_data.b0)))); elseif strcmp(data_request_eff,'rhotor')
gdat_data.dim = phi_tor.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes gdat_data.data = sqrt(phi_tor.data./pi./(ones(size(phi_tor.data,1),1)*reshape(gdat_data.b0,1,numel(gdat_data.b0))));
gdat_data.x=gdat_data.dim{1}; gdat_data.dim = phi_tor.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes
gdat_data.units = 'm'; gdat_data.x=gdat_data.dim{1};
gdat_data.dimunits = phi_tor.dimunits; gdat_data.units = 'm';
gdat_data.data_fullpath=['rhotor=sqrt(Phi/pi/B0) from phi_tor']; gdat_data.dimunits = phi_tor.dimunits;
gdat_data.data_fullpath=['rhotor=sqrt(Phi/pi/B0) from phi_tor'];
else
disp(['data_request_eff = ' data_request_eff ' not defined within rhotor block in gdat_tcv.m']);
return
end
else else
disp(['data_request_eff = ' data_request_eff ' not defined within rhotor block in gdat_tcv.m']); gdat_data.b0 = [];
return gdat_data.rhotor_edge = [];
end end
end end
...@@ -1765,9 +1774,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ...@@ -1765,9 +1774,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
psitbxput(psitbxput_version,shot); psitbxput(psitbxput_version,shot);
ishot = mdsopen(shot); ishot = mdsopen(shot);
tracetdi=tdi(nodenameeff); tracetdi=tdi(nodenameeff);
if isempty(tracetdi.data) || isempty(tracetdi.dim) end
return if isempty(tracetdi.data) || isempty(tracetdi.dim) || ischar(tracetdi.data)
end return
end end
gdat_data.units = tracetdi.units; gdat_data.units = tracetdi.units;
if strcmp(data_request_eff,'volume') if strcmp(data_request_eff,'volume')
...@@ -1788,7 +1797,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') ...@@ -1788,7 +1797,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
else else
gdat_data.data = tracetdi.data; gdat_data.data = tracetdi.data;
gdat_data.dim = tracetdi.dim; gdat_data.dim = tracetdi.dim;
gdat_data.x = gdat_data.dim{1}; if length(gdat_data.dim)>0; gdat_data.x = gdat_data.dim{1}; end
gdat_data.dimunits = tracetdi.dimunits; gdat_data.dimunits = tracetdi.dimunits;
% to always have field volume_edge % to always have field volume_edge
gdat_data.volume_edge = gdat_data.data(end,:); gdat_data.volume_edge = gdat_data.data(end,:);
......
...@@ -37,7 +37,6 @@ psi_axis = gdat(gdat_data.shot,params_eff); ...@@ -37,7 +37,6 @@ psi_axis = gdat(gdat_data.shot,params_eff);
params_eff.data_request='psi_edge'; params_eff.data_request='psi_edge';
psi_edge = gdat(gdat_data.shot,params_eff); psi_edge = gdat(gdat_data.shot,params_eff);
psi0_edge = interpos(63,psi_axis.t,psi_axis.data-psi_edge.data,gdat_data.t,-0.01); psi0_edge = interpos(63,psi_axis.t,psi_axis.data-psi_edge.data,gdat_data.t,-0.01);
if (nbdim_x == 1) if (nbdim_x == 1)
gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2*reshape(psi0_edge,1,length(psi0_edge)); gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2*reshape(psi0_edge,1,length(psi0_edge));
elseif (nbdim_x == 2) elseif (nbdim_x == 2)
...@@ -49,6 +48,11 @@ end ...@@ -49,6 +48,11 @@ end
gdat_data.grids_1d.rhotornorm = NaN*ones(size(gdat_data.data)); gdat_data.grids_1d.rhotornorm = NaN*ones(size(gdat_data.data));
gdat_data.grids_1d.rhovolnorm = NaN*ones(size(gdat_data.data)); gdat_data.grids_1d.rhovolnorm = NaN*ones(size(gdat_data.data));
if (isempty(rhotor_norm.x) ||isempty(rhotor_norm.t) || isempty(rhotor_norm.data)) ...
&& (isempty(rhovol.x) ||isempty(rhovol.t) || isempty(rhovol.data))
return
end
it_rt = iround_os(rhotor_norm.t,gdat_data.t); it_rt = iround_os(rhotor_norm.t,gdat_data.t);
it_vol = iround_os(rhovol.t,gdat_data.t); it_vol = iround_os(rhovol.t,gdat_data.t);
for it=1:length(gdat_data.t) for it=1:length(gdat_data.t)
...@@ -62,12 +66,14 @@ for it=1:length(gdat_data.t) ...@@ -62,12 +66,14 @@ for it=1:length(gdat_data.t)
end end
if (nbdim_x == 1) if (nbdim_x == 1)
if length(ii)==length(gdat_data.grids_1d.rhopolnorm) if length(ii)==length(gdat_data.grids_1d.rhopolnorm)
if ndim_x_rhotor==1 if ~isempty(rhotor_norm.x) && ~isempty(rhotor_norm.data)
gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm); if ndim_x_rhotor==1
else gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm);
gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm); else
gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm);
end
end end
if ~isempty(rhovol.x) if ~isempty(rhovol.x) && ~isempty(rhovol.data)
if ndim_x_rhovol==1 if ndim_x_rhovol==1
gdat_data.grids_1d.rhovolnorm(:,it)=interpos(-13,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm); gdat_data.grids_1d.rhovolnorm(:,it)=interpos(-13,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm);
else else
...@@ -77,12 +83,14 @@ for it=1:length(gdat_data.t) ...@@ -77,12 +83,14 @@ for it=1:length(gdat_data.t)
end end
else else
if length(ii)==size(gdat_data.grids_1d.rhopolnorm,1) if length(ii)==size(gdat_data.grids_1d.rhopolnorm,1)
if ndim_x_rhotor==1 if ~isempty(rhotor_norm.x) && ~isempty(rhotor_norm.data)
gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it)); if ndim_x_rhotor==1
else gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it));
gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it)); else
gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it));
end
end end
if ~isempty(rhovol.x) if ~isempty(rhovol.x) && ~isempty(rhovol.data)
if ndim_x_rhovol==1 if ndim_x_rhovol==1
gdat_data.grids_1d.rhovolnorm(:,it)=interpos(-13,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(:,it)); gdat_data.grids_1d.rhovolnorm(:,it)=interpos(-13,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(:,it));
else else
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment