From 51d95d3b513a64af3f29fee86aeba4eb5f2b907c Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Wed, 4 Nov 2015 17:00:31 +0000
Subject: [PATCH] start adding other cases nete_rho, cxrs, etc
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@5206 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/AUG/aug_requests_mapping.m | 2 +
crpptbx/AUG/gdat_aug.m | 209 +++++++++++++++++++++++++----
crpptbx/AUG/geteqdskAUG.m | 2 +-
crpptbx/TCV/gdat_tcv.m | 6 +-
crpptbx/gdat.m | 5 +-
5 files changed, 188 insertions(+), 36 deletions(-)
diff --git a/crpptbx/AUG/aug_requests_mapping.m b/crpptbx/AUG/aug_requests_mapping.m
index a8c81f2f..28f2a130 100644
--- a/crpptbx/AUG/aug_requests_mapping.m
+++ b/crpptbx/AUG/aug_requests_mapping.m
@@ -224,9 +224,11 @@ switch lower(data_request)
mapping.gdat_timedim = 2;
mapping.method = 'switchcase';
case 'te'
+ mapping.timedim = 2;
mapping.label = 'Te';
mapping.method = 'switchcase';
case 'te_rho'
+ mapping.timedim = 2;
mapping.label = 'Te';
mapping.method = 'switchcase';
case 'ti'
diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m
index c9ab6b69..54fe15cf 100644
--- a/crpptbx/AUG/gdat_aug.m
+++ b/crpptbx/AUG/gdat_aug.m
@@ -387,7 +387,13 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
exp_name_eff = gdat_data.gdat_params.exp_name;
diag_name = 'CEZ';
% R, Z positions of measurements
- [r_time]=sf2ab(diag_name,shot,'R_time','-exp',exp_name_eff);
+ try
+ [r_time]=sf2ab(diag_name,shot,'R_time','-exp',exp_name_eff);
+ catch ME_R_time
+ % assume no shotfile
+ getReport(ME_R_time,'basic');
+ return
+ end
gdat_data.r = r_time.value{1};
inotok=find(gdat_data.r<=0);
gdat_data.r(inotok) = NaN;
@@ -423,9 +429,11 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
[aerr,e]=rdaAUG_eff(shot,diag_name,'err_Ti_c',exp_name_eff);
gdat_data.ti.data = a.data;
gdat_data.data = a.data;
+ gdat_data.label = [gdat_data.label '/Ti'];
if any(strcmp(fieldnames(a),'units')); gdat_data.ti.units=a.units; end
if any(strcmp(fieldnames(a),'units')); gdat_data.units=a.units; end
if any(strcmp(fieldnames(a),'name')); gdat_data.ti.name=a.name; end
+ gdat_data.ti.label = 'Ti_c';
gdat_data.ti.error_bar = aerr.data;
gdat_data.error_bar = aerr.data;
gdat_data.data_fullpath=[exp_name_eff '/' diag_name '/' a.name ';vrot, Ti in data, {r;z} in .x'];
@@ -443,7 +451,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
inb_chord_cxrs=size(gdat_data.r,1);
inb_time_cxrs=size(gdat_data.r,2);
psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
- rhopsinorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
+ rhopolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
@@ -459,19 +467,20 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
zout=gdat_data.z(iok,it_cxrs_inequil);
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil));
- rhopsinorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
+ rhopolnorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
for it_cx=1:length(it_cxrs_inequil)
- rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
- rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
+ rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
+ rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
end
gdat_data.psi = psi_out;
- gdat_data.rhopsinorm = rhopsinorm_out;
+ gdat_data.rhopolnorm = rhopolnorm_out;
gdat_data.rhotornorm = rhotornorm_out;
gdat_data.rhovolnorm = rhovolnorm_out;
end
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ece', 'eced', 'ece_rho', 'eced_rho'}
nth_points = 13;
if isfield(gdat_data.gdat_params,'nth_points') && ~isempty(gdat_data.gdat_params.nth_points)
@@ -573,7 +582,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
inb_chord_ece=size(gdat_data.rz.r,1);
inb_time_ece=size(gdat_data.rz.r,2);
psi_out = NaN*ones(inb_chord_ece,inb_time_ece);
- rhopsinorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
+ rhopolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
rhotornorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
rhovolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
@@ -595,21 +604,22 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
for ij=1:length(iok)
psi_out(ieff(ij),it_ece_inequil(jeff(ij))) = psi_at_routzout(ij);
end
- rhopsinorm_out(:,it_ece_inequil) = sqrt(abs((psi_out(:,it_ece_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))));
+ rhopolnorm_out(:,it_ece_inequil) = sqrt(abs((psi_out(:,it_ece_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))));
for it_cx=1:length(it_ece_inequil)
- rhotornorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
- rhovolnorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
+ rhotornorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
+ rhovolnorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
end
end
gdat_data.rhos.psi = psi_out;
- gdat_data.rhos.rhopsinorm = rhopsinorm_out;
+ gdat_data.rhos.rhopolnorm = rhopolnorm_out;
gdat_data.rhos.rhotornorm = rhotornorm_out;
gdat_data.rhos.rhovolnorm = rhovolnorm_out;
gdat_data.rhos.t = gdat_data.rz.t;
end
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'eqdsk'}
%
time_eqdsks=1.; % default time
@@ -684,6 +694,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.request_description = ['data=psi, x=(R,Z), eqdsk contains eqdsk structure with which ' ...
'plot_eqdsk, write_eqdsk, read_eqdsk can be used'];
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'equil'}
exp_name_eff = gdat_data.gdat_params.exp_name;
gdat_data.gdat_params.equil = upper(gdat_data.gdat_params.equil);
@@ -878,9 +889,10 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.cocos = 17; % should check FPP
gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
- gdat_data.dimunits{1} = '';
- gdat_data.dimunits{2} = 's';
+ gdat_data.dimunits{1} = 'rhopolnorm';
+ gdat_data.dimunits{2} = 'time [s]';
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ne','te'}
exp_name_eff = gdat_data.gdat_params.exp_name;
% ne or Te from Thomson data on raw z mesh vs (z,t)
@@ -946,9 +958,10 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
end
gdat_data.x=gdat_data.dim{1};
- case {'ne_rho', 'te_rho'}
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'ne_rho', 'te_rho', 'nete_rho'}
params_eff = gdat_data.gdat_params;
- params_eff.data_request=data_request_eff(1:2);
+ params_eff.data_request=data_request_eff(1:2); % start with ne if nete_rho
% get raw data
[gdat_data,params_kin,error_status]=gdat_aug(shot,params_eff);
if error_status>0
@@ -970,7 +983,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
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);
- rhopsinorm_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
@@ -978,7 +991,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
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);
- rhopsinorm_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)]
@@ -993,12 +1006,12 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
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));
- rhopsinorm_out_core(:,it_core_inequil) = sqrt((psi_out_core(:,it_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
+ 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),rhopsinorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
+ 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),rhopsinorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
+ interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
% edge
@@ -1008,31 +1021,169 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
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));
- rhopsinorm_out_edge(:,it_edge_inequil) = sqrt((psi_out_edge(:,it_edge_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
+ 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),rhopsinorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
+ 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),rhopsinorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
+ interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]);
end
end
end
gdat_data.(lower(node_child_nameeff)).psi = psi_out_core;
- gdat_data.(lower(node_child_nameeff)).rhopsinorm = rhopsinorm_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)).rhopsinorm = rhopsinorm_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 rhopsinorm for dim{1} by default
- gdat_data.x = gdat_data.(lower(node_child_nameeff)).rhopsinorm;
+ % put values of rhopolnorm for dim{1} by default
+ gdat_data.x = gdat_data.(lower(node_child_nameeff)).rhopolnorm;
iaaa=iroundos(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)).rhopsinorm(:,iaaa);
+ 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} = 'rhopsinorm';
+ gdat_data.dimunits{1} = 'rhopolnorm';
gdat_data.data_fullpath = [gdat_data.data_fullpath ' projected on equilibrium ' gdat_data.gdat_params.equil];
+
+ % if nete_rho, copy data as .ne, get .te and put pe=e ne Te in data:
+ if strcmp(data_request_eff(1:4),'nete')
+ % note, now has ne.data_raw for density without fir_to_thomson_ratio
+ gdat_data.ne.data = gdat_data.data;
+ gdat_data.ne.error_bar = gdat_data.error_bar;
+ gdat_data.ne.units = 'm^{-3}';
+ gdat_data.ne.ne_core = gdat_data.ne_core;
+ gdat_data.ne.ne_edge = gdat_data.ne_edge;
+ gdat_data = rmfield(gdat_data,{'ne_core','ne_edge'});
+ %
+ params_eff.data_request=data_request_eff(3:4);
+ [gdat_data_te,params_kin,error_status]=gdat_aug(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.te_core = gdat_data_te.te_core;
+ gdat_data.te.te_edge = gdat_data_te.te_edge;
+ gdat_data.te.error_bar = gdat_data_te.te.error_bar;
+ gdat_data.te.te_core.psi = gdat_data.ne.ne_core.psi;
+ gdat_data.te.te_core.rhopolnorm = gdat_data.ne.ne_core.rhopolnorm;
+ gdat_data.te.te_core.rhotornorm = gdat_data.ne.ne_core.rhotornorm;
+ gdat_data.te.te_core.rhovolnorm = gdat_data.ne.ne_core.rhovolnorm;
+ gdat_data.te.te_edge.psi = gdat_data.ne.ne_edge.psi;
+ gdat_data.te.te_edge.rhopolnorm = gdat_data.ne.ne_edge.rhopolnorm;
+ gdat_data.te.te_edge.rhotornorm = gdat_data.ne.ne_edge.rhotornorm;
+ gdat_data.te.te_edge.rhovolnorm = gdat_data.ne.ne_edge.rhovolnorm;
+ end
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'q_rho'}
+ exp_name_eff = gdat_data.gdat_params.exp_name;
+ gdat_data.gdat_params.equil = upper(gdat_data.gdat_params.equil);
+ DIAG = gdat_data.gdat_params.equil; % 'EQI' by default at this stage, should become EQH?
+ M_Rmesh_par = sf2par(DIAG,shot,'M','PARMV');
+ M_Rmesh = M_Rmesh_par.value + 1; % nb of points
+ N_Zmesh_par = sf2par(DIAG,shot,'N','PARMV');
+ N_Zmesh = N_Zmesh_par.value + 1; % nb of points
+ NTIME_par = sf2par(DIAG,shot,'NTIME','PARMV');
+ NTIME = NTIME_par.value; % nb of points
+ Lpf_par = rdaAUG_eff(shot,DIAG,'Lpf',exp_name_eff);
+ % since June, nb of time points in EQ results is not consistent with NTIME and time
+ % It seems the first NTIME points are correct, so use this explicitely
+ NTIME_Lpf = length(Lpf_par.value);
+ if (NTIME < NTIME_Lpf)
+ if gdat_params.nverbose>=3; disp('WARNING: nb of times points smaller then equil results, use first NTIME points'); end
+ elseif (NTIME > NTIME_Lpf)
+ if gdat_params.nverbose >= 1
+ disp('ERROR: nb of times points LARGER then equil results')
+ disp('this is unexpected, so stop there and ask Olivier.Sauter@epfl.ch')
+ end
+ return
+ end
+ Lpf_tot = Lpf_par.value(1:NTIME); % nb of points: 100000*nb_SOL points + nb_core
+ Lpf_SOL = fix(Lpf_tot/100000);
+ Lpf1_t = mod(Lpf_tot,100000)+1; % nb of Lpf points
+ % since Lpf depends on time, need to load all first and then loop over time for easier mapping
+ [qpsi,e]=rdaAUG_eff(shot,DIAG,'Qpsi',exp_name_eff);
+ ndimrho = size(qpsi.data,2);
+ if ndimrho==NTIME_Lpf
+ % data seems to be transposed
+ ndimrho = size(qpsi.data,1);
+ itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
+ else
+ itotransposeback = 0;
+ end
+ qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback);
+ ijnan=find(isnan(qpsi.value));
+ qpsi.value(ijnan)=0;
+ [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff);
+ psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
+ [phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',exp_name_eff);
+ phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
+ [Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',exp_name_eff);
+ Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
+ % seems "LCFS" q-value is far too large, limit to some max (when diverted)
+ max_qValue = 30.; % Note could just put a NaN on LCFS value since ill-defined when diverted
+ for it=1:NTIME
+ Lpf1 = Lpf1_t(it);
+ % get x values
+ psi_it=psi_tree.value(it,Lpf1:-1:1)';
+ gdat_data.psi_axis(it)= psi_it(1);
+ gdat_data.psi_lcfs(it)= psi_it(end);
+ gdat_data.rhopolnorm(:,it) = sqrt(abs((psi_it-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it))));
+ if strcmp(DIAG,'EQR');
+ % q value has only a few values and from center to edge, assume they are from central rhopol values on
+ % But they are every other point starting from 3rd
+ ijk=find(gdat_data.qvalue(:,it)~=0);
+ if length(ijk)>2
+ % now shots have non-zero axis values in eqr
+ rhoeff=gdat_data.rhopolnorm(ijk,it);
+ qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above
+ if ijk(1)>1
+ rhoeff = [0.; rhoeff];
+ qeff = [qeff(1) ;qeff];
+ end
+ ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it)));
+ qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
+ qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]);
+ else
+ qfit = zeros(size(gdat_data.rhopolnorm(:,it)));
+ end
+ gdat_data.qvalue(:,it) = qfit;
+ end
+ % get rhotor values
+ phi_it = phi_tree.value(it,Lpf1:-1:1)';
+ gdat_data.rhotornorm(:,it) = sqrt(abs(phi_it ./ phi_it(end)));
+ % get rhovol values
+ vol_it=Vol.value(it,2*Lpf1-1:-2:1)';
+ gdat_data.rhovolnorm(:,it) = sqrt(abs(vol_it ./ vol_it(end)));
+ % Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
+ % change it to (radial,time) and use only Lpf+1 points up to LCFS
+ ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially
+ % max value 1e6, change to 2*extrapolation (was max_qValue before)
+ q_edge=interpos(gdat_data.rhotornorm(1:end-1,it),qpsi.value(it,Lpf1:-1:2),1,-0.1);
+ gdat_data.qvalue(:,it) = qpsi.value(it,Lpf1:-1:1)';
+ if abs(gdat_data.qvalue(end,it)) > 1e3
+ % assume diverted
+ 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);
+% $$$ else
+% $$$ gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
+% $$$ end
+ end
+ gdat_data.x = gdat_data.rhopolnorm;
+ % get time values
+ [equil_time,e]=rdaAUG_eff(shot,DIAG,'time',exp_name_eff);
+ gdat_data.t = equil_time.value(1:NTIME);
+ gdat_data.data = gdat_data.qvalue; % put q in data
+ gdat_data.units=[]; % not applicable
+ gdat_data.data_fullpath = [DIAG ' from expname: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
+ gdat_data.cocos = 17; % should check FPP
+ gdat_data.dim{1} = gdat_data.x;
+ gdat_data.dim{2} = gdat_data.t;
+ gdat_data.dimunits{1} = 'rhopolnorm';
+ gdat_data.dimunits{2} = 'time [s]';
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'sxr'}
% sxr from B by default or else if 'source','else' is provided
diff --git a/crpptbx/AUG/geteqdskAUG.m b/crpptbx/AUG/geteqdskAUG.m
index fc6a8a3c..7a661a18 100644
--- a/crpptbx/AUG/geteqdskAUG.m
+++ b/crpptbx/AUG/geteqdskAUG.m
@@ -87,7 +87,7 @@ eqdsk.zaxis = zmag.data(itrmag) - eqdsk.zshift;
% get plasma boundary
figure
contour(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',100)
-hold
+hold on
psiedge_eff = 0.001*eqdsk.psiaxis + 0.999*eqdsk.psiedge;
[hh1 hh2]=contour(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',[psiedge_eff psiedge_eff],'k');
axis equal
diff --git a/crpptbx/TCV/gdat_tcv.m b/crpptbx/TCV/gdat_tcv.m
index b4fc6be1..879b38a4 100644
--- a/crpptbx/TCV/gdat_tcv.m
+++ b/crpptbx/TCV/gdat_tcv.m
@@ -1335,7 +1335,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
% invert index of time and channel (rho)
gdat_data.data = mpx.(gdat_data.gdat_params.camera(1:3)).signal.data';
gdat_data.t = mpx.(gdat_data.gdat_params.camera(1:3)).signal.dim{1};
- gdat_data.dim{1} = interp1(mpx.top.rho.time,mpx.top.rho.rhopsi,gdat_data.t);
+ gdat_data.dim{1} = interp1(mpx.top.rho.time,mpx.top.rho.rhopol,gdat_data.t);
gdat_data.dim{2} = mpx.(gdat_data.gdat_params.camera(1:3)).signal.dim{1};
gdat_data.x = gdat_data.dim{1};
gdat_data.(gdat_data.gdat_params.camera).data = gdat_data.data;
@@ -1346,14 +1346,14 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
else
gdat_data.data = mpx.signal.data';
gdat_data.t = mpx.signal.dim{1};
- gdat_data.dim{1} = interp1(mpx.top.rho.time,mpx.top.rho.rhopsi,gdat_data.t);
+ gdat_data.dim{1} = interp1(mpx.top.rho.time,mpx.top.rho.rhopol,gdat_data.t);
gdat_data.dim{2} = mpx.(gdat_data.gdat_params.camera(1:3)).signal.dim{1};
%
gdat_data.top.data = mpx.top.signal.data';
gdat_data.top.x = gdat_data.dim{1};
gdat_data.top.channel = mpx.top.signal.dim{2};
gdat_data.bottom.data = mpx.bot.signal.data;
- gdat_data.bottom.x = interp1(mpx.bot.rho.time,mpx.bot.rho.rhopsi,gdat_data.t);
+ gdat_data.bottom.x = interp1(mpx.bot.rho.time,mpx.bot.rho.rhopol,gdat_data.t);
gdat_data.bottom.channel = mpx.bottom.signal.dim{2};
gdat_data.(gdat_data.gdat_params.camera).channel = mpx.(gdat_data.gdat_params.camera(1:3)).signal.dim{2};
gdat_data.data_fullpath = ['MPX for ' gdat_data.gdat_params.camera ' camera in .data, "rho" in .x between [-1,1]' ...
diff --git a/crpptbx/gdat.m b/crpptbx/gdat.m
index 7291d3ce..7566bbe4 100644
--- a/crpptbx/gdat.m
+++ b/crpptbx/gdat.m
@@ -170,7 +170,7 @@ catch ME_gdat
if ~exist('gdat_params'); gdat_params.plot = []; end
if ~exist('error_status'); error_status = 998; end
if exist('ME_gdat')
- rethrow(ME_gdat)
+ getReport(ME_gdat)
end
return
end
@@ -184,7 +184,6 @@ if gdat_data.gdat_params.doplot
catch ME_gdat_plot
hhh = [];
end
-5
if nargout >=3
if length(varargout)==0 || isempty(varargout{1})
varargout{1} = hhh;
@@ -195,5 +194,5 @@ if gdat_data.gdat_params.doplot
end
if exist('ME_gdat_plot')
- rethrow(ME_gdat_plot)
+ getReport(ME_gdat_plot)
end
--
GitLab