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

add several D3D signals, including gdat(shot,'eqdsk') and sxr, ece 

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@7299 d63d8f72-b253-0410-a779-e742ad2e26cf
parent 5fe3f50b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
crpptbx/D3D/d3d_requests_mapping.m
+ 23
− 18
View file @ cde58f23
...@@ -52,7 +52,7 @@ switch lower(data_request) ...@@ -52,7 +52,7 @@ switch lower(data_request)
mapping.timedim = 1; mapping.timedim = 1;
mapping.label = 'B_0'; mapping.label = 'B_0';
mapping.method = 'signal'; mapping.method = 'signal';
mapping.expression = [{'FPC'},{'BTF'}]; mapping.expression = [{'EFIT01'},{'\bcentr'}];
case 'beta' case 'beta'
mapping.timedim = 1; mapping.timedim = 1;
mapping.label = '\beta'; mapping.label = '\beta';
...@@ -106,23 +106,23 @@ switch lower(data_request) ...@@ -106,23 +106,23 @@ switch lower(data_request)
mapping.label = 'delta\_top'; mapping.label = 'delta\_top';
mapping.timedim = 1; mapping.timedim = 1;
mapping.method = 'signal'; mapping.method = 'signal';
mapping.expression = [{'FPG'},{'delRoben'}]; mapping.expression = [{'OT01'},{'\tritop'}];
case 'delta_bottom' case 'delta_bottom'
mapping.label = 'delta\_bottom'; mapping.label = 'delta\_bottom';
mapping.timedim = 1; mapping.timedim = 1;
mapping.method = 'signal'; mapping.method = 'signal';
mapping.expression = [{'FPG'},{'delRuntn'}]; mapping.expression = [{'OT01'},{'\tribot'}];
case {'ece', 'eced', 'ece_rho', 'eced_rho'} case {'ece', 'ece_rho'}
mapping.timedim = 2; mapping.timedim = 1;
mapping.method = 'switchcase'; mapping.method = 'switchcase';
mapping.expression = ''; mapping.expression = '';
case 'eqdsk' case 'eqdsk'
mapping.timedim = 2; mapping.timedim = 2;
mapping.method = 'switchcase'; % could use function make_eqdsk directly? mapping.method = 'switchcase';
mapping.expression = ''; mapping.expression = '';
case 'equil' case 'equil'
mapping.gdat_timedim = 2; mapping.gdat_timedim = 2;
mapping.method = 'switchcase'; % could use function make_eqdsk directly? mapping.method = 'switchcase';
mapping.expression = ''; mapping.expression = '';
case 'halpha' case 'halpha'
mapping.timedim = 1; mapping.timedim = 1;
...@@ -148,12 +148,12 @@ switch lower(data_request) ...@@ -148,12 +148,12 @@ switch lower(data_request)
mapping.timedim = 1; mapping.timedim = 1;
mapping.label = 'Plasma current'; mapping.label = 'Plasma current';
mapping.method = 'signal'; mapping.method = 'signal';
mapping.expression = [{'MAG'},{'Ipa'}]; mapping.expression = [{'efit01'},{'\cpasma'}];
case 'kappa' case 'kappa'
mapping.timedim = 1; mapping.timedim = 1;
mapping.label = '\kappa'; mapping.label = '\kappa';
mapping.method = 'signal'; mapping.method = 'signal';
mapping.expression = [{'FPG'},{'k'}]; mapping.expression = [{'efit01'},{'\kappa'}];
case 'kappa_top' case 'kappa_top'
mapping.timedim = 1; mapping.timedim = 1;
mapping.label = '\kappa^{top}'; mapping.label = '\kappa^{top}';
...@@ -168,7 +168,7 @@ switch lower(data_request) ...@@ -168,7 +168,7 @@ switch lower(data_request)
mapping.timedim = 1; mapping.timedim = 1;
mapping.label = 'l_i'; mapping.label = 'l_i';
mapping.method = 'signal'; mapping.method = 'signal';
mapping.expression = [{'FPG'},{'li'}]; mapping.expression = [{'efit01'},{'\li3'}];
case 'mhd' case 'mhd'
mapping.timedim = 1; mapping.timedim = 1;
mapping.label = 'Odd and Even n'; mapping.label = 'Odd and Even n';
...@@ -188,6 +188,11 @@ switch lower(data_request) ...@@ -188,6 +188,11 @@ switch lower(data_request)
'gdat_tmp.n_even.amp_t=gdat_tmp2.t;' ... 'gdat_tmp.n_even.amp_t=gdat_tmp2.t;' ...
'gdat_tmp.full_path=''MOD/Odd in data and .n_odd; .n_even'';' ... 'gdat_tmp.full_path=''MOD/Odd in data and .n_odd; .n_even'';' ...
'gdat_tmp.gdat_request=''mhd'';gdat_tmp.gdat_params.data_request=gdat_tmp.gdat_request;']; 'gdat_tmp.gdat_request=''mhd'';gdat_tmp.gdat_params.data_request=gdat_tmp.gdat_request;'];
case 'nbi'
mapping.timedim = 1;
mapping.label = 'Pnbi';
mapping.method = 'signal';
mapping.expression = [{'nb'},{'pinj'}];
case 'ne' case 'ne'
mapping.timedim = 2; mapping.timedim = 2;
mapping.method = 'switchcase'; mapping.method = 'switchcase';
...@@ -198,14 +203,9 @@ switch lower(data_request) ...@@ -198,14 +203,9 @@ switch lower(data_request)
mapping.expression = [{'DCN'},{'H-1'}]; mapping.expression = [{'DCN'},{'H-1'}];
case 'nel' case 'nel'
mapping.timedim = 1; mapping.timedim = 1;
mapping.label = 'line-averaged el. density'; mapping.label = 'NEBAR\_R0';
mapping.expression = [{'FPG'},{'lenH-1'}]; mapping.method = 'signal';
mapping.method = 'expression'; mapping.expression = [{'efit01'},{'\NEBAR_R0'}];
mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''neint'';' ...
'gdat_tmp=gdat_d3d(shot,params_eff);params_eff.data_request=[{''FPG''},{''lenH-1''}];' ...
'gdat_tmp2=gdat_d3d(shot,params_eff);ij=find(gdat_tmp2.data==0);gdat_tmp2.data(ij)=NaN;' ...
'tmp_data=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ...
'gdat_tmp.data = gdat_tmp.data./(tmp_data+1e-5);'];
case 'ne_rho' case 'ne_rho'
mapping.timedim = 2; mapping.timedim = 2;
mapping.label = 'ne'; mapping.label = 'ne';
...@@ -254,6 +254,11 @@ switch lower(data_request) ...@@ -254,6 +254,11 @@ switch lower(data_request)
mapping.gdat_timedim = 2; mapping.gdat_timedim = 2;
mapping.label = 'q'; mapping.label = 'q';
mapping.method = 'switchcase'; mapping.method = 'switchcase';
case 'r0'
mapping.timedim = 1;
mapping.label = 'R_0';
mapping.method = 'signal';
mapping.expression = [{'EFIT01'},{'\rzero'}];
case 'rgeom' case 'rgeom'
mapping.label = 'Rgeom'; mapping.label = 'Rgeom';
mapping.timedim = 1; mapping.timedim = 1;
......
crpptbx/D3D/gdat_d3d.m
+ 145
− 470
View file @ cde58f23
...@@ -198,9 +198,7 @@ if (nargin>=ivarargin_first_char) ...@@ -198,9 +198,7 @@ if (nargin>=ivarargin_first_char)
end end
data_request_eff = gdat_params.data_request; % in case was defined in pairs data_request_eff = gdat_params.data_request; % in case was defined in pairs
% default equil: % default equil:
if ~isfield(gdat_params,'equil'); gdat_params.equil = 'EQI'; end if ~isfield(gdat_params,'equil'); gdat_params.equil = 'EFIT01'; end
if isfield(gdat_params,'source') && (any(strcmp(lower(gdat_params.source),'ide')) || any(strcmp(lower(gdat_params.source),'idg')) ...
|| any(strcmp(lower(gdat_params.source),'ida'))); gdat_params.equil = 'IDE'; end
% if it is a request_keyword can obtain description: % if it is a request_keyword can obtain description:
if ischar(data_request_eff) || length(data_request_eff)==1 if ischar(data_request_eff) || length(data_request_eff)==1
...@@ -271,7 +269,9 @@ if strcmp(mapping_for_d3d.method,'signal') ...@@ -271,7 +269,9 @@ if strcmp(mapping_for_d3d.method,'signal')
if gdat_params.nverbose>=3; disp(['error after get_signal_d3d in signal with data_request_eff= ' data_request_eff]); end if gdat_params.nverbose>=3; disp(['error after get_signal_d3d in signal with data_request_eff= ' data_request_eff]); end
return return
end end
if iscell(gdat_data.label) && length(gdat_data.label)==2;
gdat_data.label = {[gdat_data.label{1} ' ' gdat_data.label{2}]};
end
gdat_data.data_fullpath=mapping_for_d3d.expression; gdat_data.data_fullpath=mapping_for_d3d.expression;
% end of method "signal" % end of method "signal"
...@@ -533,7 +533,7 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') ...@@ -533,7 +533,7 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
end end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'ece', 'eced', 'ece_rho', 'eced_rho'} case {'ece', 'ece_rho'}
nth_points = 13; nth_points = 13;
if isfield(gdat_data.gdat_params,'nth_points') && ~isempty(gdat_data.gdat_params.nth_points) if isfield(gdat_data.gdat_params,'nth_points') && ~isempty(gdat_data.gdat_params.nth_points)
nth_points = gdat_data.gdat_params.nth_points; nth_points = gdat_data.gdat_params.nth_points;
...@@ -544,6 +544,7 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') ...@@ -544,6 +544,7 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels) if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels)
channels = gdat_data.gdat_params.channels; channels = gdat_data.gdat_params.channels;
end end
diag_name = 'ece';
if nth_points>=10 if nth_points>=10
match_rz_to_time = 1; match_rz_to_time = 1;
else else
...@@ -561,126 +562,66 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') ...@@ -561,126 +562,66 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
gdat_data.gdat_params.time_interval = time_interval; gdat_data.gdat_params.time_interval = time_interval;
end end
% %
if isfield(gdat_data.gdat_params,'diag_name') && ~isempty(gdat_data.gdat_params.diag_name)
diag_name = gdat_data.gdat_params.diag_name;
if strcmp(upper(diag_name),'RMD'); gdat_data.gdat_params.exp_name = 'ECED'; end
else
diag_name = 'CEC';
gdat_data.gdat_params.diag_name = diag_name;
end
exp_name_eff = gdat_data.gdat_params.exp_name;
if strcmp(data_request_eff,'eced')
exp_name_eff = 'ECED';
gdat_data.gdat_params.exp_name = exp_name_eff;
end
[a,e]=rdaD3D_eff(shot,diag_name,'Trad-A',exp_name_eff,time_interval);
% [a,e]=rdaD3D_eff(shot,diag_name,'Trad-A',exp_name_eff);
inb_chord = size(a.data,1);
if channels(1)<=0 if channels(1)<=0
channels = [1:inb_chord]; channels = [1:40];
end end
chanelles=sort(channels);
gdat_data.dim{1} = channels; gdat_data.dim{1} = channels;
gdat_data.gdat_params.channels = channels; gdat_data.gdat_params.channels = channels;
if nth_points>1 for i=channels
gdat_data.data = a.data(channels,1:nth_points:end); a = gdat_d3d(shot,{'ece',['\tece' num2str(i,'%.2d')]});
gdat_data.dim{2} = a.t(1:nth_points:end); if nth_points>1
else gdat_data.data(i,:) = a.data(1:nth_points:end);
gdat_data.data = a.data(channels,:); gdat_data.dim{2} = a.t(1:nth_points:end);
gdat_data.dim{2} = a.t; else
gdat_data.data(i,:) = a.data;
gdat_data.dim{2} = a.t;
end
end end
gdat_data.x = gdat_data.dim{1}; gdat_data.x = gdat_data.dim{1};
gdat_data.t = gdat_data.dim{2}; gdat_data.t = gdat_data.dim{2};
gdat_data.dimunits=[{'channels'} ; {'time [s]'}]; gdat_data.dimunits=[{'channels'} ; {'time [s]'}];
if any(strcmp(fieldnames(a),'units')); gdat_data.units=a.units; end if any(strcmp(fieldnames(a),'units')); gdat_data.units=a.units; end
try % $$$ if match_rz_to_time
[aR,e]=rdaD3D_eff(shot,diag_name,'R-A',exp_name_eff,time_interval); % $$$ % interpolate R structure on ece data time array, to ease plot vs R
catch % $$$ for i=1:length(channels)
end % $$$ radius.data(i,:) = interp1(aR.t,aR.data(channels(i),:),gdat_data.t);
try % $$$ zheight.data(i,:) = interp1(aZ.t,aZ.data(channels(i),:),gdat_data.t);
[aZ,e]=rdaD3D_eff(shot,diag_name,'z-A',exp_name_eff,time_interval); % $$$ end
catch % $$$ radiuszheight.t=gdat_data.t;
disp(['problem with getting z-A in ' diag_name]) % $$$ else
end % $$$ radius.data = aR.data(channels,:);
if match_rz_to_time % $$$ radiuszheight.t=aR.t;
% interpolate R structure on ece data time array, to ease plot vs R % $$$ zheight.data = aZ.data(channels,:);
for i=1:length(channels) % $$$ end
radius.data(i,:) = interp1(aR.t,aR.data(channels(i),:),gdat_data.t);
zheight.data(i,:) = interp1(aZ.t,aZ.data(channels(i),:),gdat_data.t);
end
radiuszheight.t=gdat_data.t;
else
radius.data = aR.data(channels,:);
radiuszheight.t=aR.t;
zheight.data = aZ.data(channels,:);
end
ij=find(gdat_data.data<=0); ij=find(gdat_data.data<=0);
gdat_data.data(ij)=NaN; gdat_data.data(ij)=NaN;
gdat_data.rz.r=radius.data; % $$$ gdat_data.rz.r=radius.data;
gdat_data.rz.z=zheight.data; % $$$ gdat_data.rz.z=zheight.data;
gdat_data.rz.t = radiuszheight.t; % $$$ gdat_data.rz.t = radiuszheight.t;
gdat_data.data_fullpath = [exp_name_eff '/' diag_name '/Trad-A with R, Z in .r,.z']; gdat_data.data_fullpath = [diag_name '\tece' num2str(channels(1)) '...' num2str(channels(end))];
if strcmp(data_request_eff,'ece_rho') || strcmp(data_request_eff,'eced_rho') % $$$ gdat_data.rhos.psi = psi_out;
params_equil = gdat_data.gdat_params; % $$$ gdat_data.rhos.rhopolnorm = rhopolnorm_out;
params_equil.data_request = 'equil'; % $$$ gdat_data.rhos.rhotornorm = rhotornorm_out;
[equil,params_equil,error_status] = gdat_d3d(shot,params_equil); % $$$ gdat_data.rhos.rhovolnorm = rhovolnorm_out;
if error_status>0 % $$$ gdat_data.rhos.t = gdat_data.rz.t;
if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end % $$$ end
return
end
gdat_data.gdat_params.equil = params_equil.equil;
gdat_data.equil = equil;
gdat_data.data_fullpath = [exp_name_eff '/' diag_name '/Trad-A with .r,.z projected on equil ' gdat_data.gdat_params.equil ' in .rhos'];
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);
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)]
time_equil=[min(gdat_data.rz.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.rz.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_ece_inequil = find(gdat_data.rz.t>time_equil(itequil) & gdat_data.rz.t<=time_equil(itequil+1));
if ~isempty(it_ece_inequil)
r_it_ece_inequil = gdat_data.rz.r(:,it_ece_inequil);
iok = find(~isnan(r_it_ece_inequil) & r_it_ece_inequil>0);
if ~isempty(iok)
rout=r_it_ece_inequil(iok);
z_it_ece_inequil = gdat_data.rz.z(:,it_ece_inequil);
zout=z_it_ece_inequil(iok);
psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
[ieff,jeff]=ind2sub(size(gdat_data.rz.r(:,it_ece_inequil)),iok);
for ij=1:length(iok)
psi_out(ieff(ij),it_ece_inequil(jeff(ij))) = psi_at_routzout(ij);
end
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),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.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'} case {'eqdsk'}
% %
time_eqdsks=1.; % default time time_eqdsks=2000.; % default time
if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time) if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time)
time_eqdsks = gdat_data.gdat_params.time; time_eqdsks = gdat_data.gdat_params.time;
else else
gdat_data.gdat_params.time = time_eqdsks; gdat_data.gdat_params.time = time_eqdsks;
disp(['"time" is expected as an option, choose default time = ' num2str(time_eqdsks)]); disp(['"time" is expected as an option, choose default time = ' num2str(time_eqdsks)]);
end end
if time_eqdsks < 10;
% assume given in seconds, d3d time is in ms
time_eqdsks = time_eqdsks * 1e3;
end
gdat_data.gdat_params.time = time_eqdsks; gdat_data.gdat_params.time = time_eqdsks;
gdat_data.t = time_eqdsks; gdat_data.t = time_eqdsks;
zshift = 0.; zshift = 0.;
...@@ -702,7 +643,7 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') ...@@ -702,7 +643,7 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
for itime=1:length(time_eqdsks) for itime=1:length(time_eqdsks)
time_eff = time_eqdsks(itime); time_eff = time_eqdsks(itime);
% use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2 % use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2
[eqdskD3D, equil_all_t, equil_t_index]=geteqdskD3D(equil,time_eff,gdat_data.gdat_params.zshift,'source',gdat_data.gdat_params.equil); [eqdskD3D, equil_all_t, equil_t_index]=get_eqdsk_d3d(equil,time_eff,gdat_data.gdat_params.zshift,'source',gdat_data.gdat_params.equil);
eqdskD3D.fnamefull = fullfile(['/tmp/' getenv('USER')],['EQDSK_' num2str(shot) 't' num2str(time_eff,'%.4f')]); eqdskD3D.fnamefull = fullfile(['/tmp/' getenv('USER')],['EQDSK_' num2str(shot) 't' num2str(time_eff,'%.4f')]);
cocos_out = equil.cocos; cocos_out = equil.cocos;
if isfield(gdat_data.gdat_params,'cocos') && ~isempty(gdat_data.gdat_params.cocos) if isfield(gdat_data.gdat_params,'cocos') && ~isempty(gdat_data.gdat_params.cocos)
...@@ -748,183 +689,70 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') ...@@ -748,183 +689,70 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'equil'} case {'equil'}
% get equil params and time array in gdat_data.t if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || ~ischar(gdat_data.gdat_params.source)
[gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data); gdat_data.gdat_params.source = 'efit03';
% since Lpf depends on time, need to load all first and then loop over time for easier mapping
[qpsi,e]=rdaD3D_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 else
itotransposeback = 0;
end end
qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback); if exist('read_mds_eq_func')<=0
ijnan=find(isnan(qpsi.value)); if gdat_params.nverbose>=3; disp(['addpath(''/m/GAtools/matlab/efit'',''end'') since needs function read_mds_eq_func']); end
qpsi.value(ijnan)=0; addpath('/m/GAtools/matlab/efit','end');
[psi_tree,e]=rdaD3D_eff(shot,DIAG,'PFL',exp_name_eff); end
psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback); [equil_all_t,neq,eq,ier]=read_mds_eq_func(shot,gdat_data.gdat_params.source,'d3d');
[phi_tree,e]=rdaD3D_eff(shot,DIAG,'TFLx',exp_name_eff); gdat_data.rhovolnorm = eq.results.geqdsk.rhovn;
phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback); gdat_data.equil_all_t = equil_all_t;
[Vol,e]=rdaD3D_eff(shot,DIAG,'Vol',exp_name_eff); gdat_data.eq_raw = eq;
Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback); % extract main parameters, can add more later
[Area,e]=rdaD3D_eff(shot,DIAG,'Area',exp_name_eff); gdat_data.phi = [];
Area=adapt_rda(Area,NTIME,2*ndimrho,itotransposeback); gdat_data.rhotornorm = [];
[Ri,e]=rdaD3D_eff(shot,DIAG,'Ri',exp_name_eff); gdat_data.t = equil_all_t.time;
Ri=adapt_rda(Ri,NTIME,M_Rmesh,itotransposeback); for it=1:length(equil_all_t.time)
[Zj,e]=rdaD3D_eff(shot,DIAG,'Zj',exp_name_eff); gdat_data.qvalue(:,it) = equil_all_t.gdata(it).qpsi;
Zj=adapt_rda(Zj,NTIME,N_Zmesh,itotransposeback); gdat_data.psi(:,it) = (equil_all_t.gdata(it).ssibry - equil_all_t.gdata(it).ssimag) ...
[PFM_tree,e]=rdaD3D_eff(shot,DIAG,'PFM','-raw','-exp',exp_name_eff); % -raw necessary for IDE .* linspace(0,1,equil_all_t.gdata(it).nw)' + equil_all_t.gdata(it).ssimag;
PFM_tree=adaptPFM_rda(PFM_tree,M_Rmesh,N_Zmesh,NTIME); gdat_data.psi_axis(it)= equil_all_t.gdata(it).ssimag;
[Pres,e]=rdaD3D_eff(shot,DIAG,'Pres',exp_name_eff); gdat_data.psi_lcfs(it)= equil_all_t.gdata(it).ssibry;
Pres=adapt_rda(Pres,NTIME,2*ndimrho,itotransposeback);
[Jpol,e]=rdaD3D_eff(shot,DIAG,'Jpol',exp_name_eff);
Jpol=adapt_rda(Jpol,NTIME,2*ndimrho,itotransposeback);
[FFP,e]=rdaD3D_eff(shot,DIAG,'FFP',exp_name_eff);
if ~isempty(FFP.value)
FFP=adapt_rda(FFP,NTIME,ndimrho,itotransposeback);
else
FFP.value=NaN*ones(NTIME,max(Lpf1_t));
end
if strcmp(DIAG,'EQI') || strcmp(DIAG,'EQH') || strcmp(DIAG,'IDE')
[Rinv,e]=rdaD3D_eff(shot,DIAG,'Rinv',exp_name_eff);
Rinv=adapt_rda(Rinv,NTIME,ndimrho,itotransposeback);
[R2inv,e]=rdaD3D_eff(shot,DIAG,'R2inv',exp_name_eff);
R2inv=adapt_rda(R2inv,NTIME,ndimrho,itotransposeback);
[Bave,e]=rdaD3D_eff(shot,DIAG,'Bave',exp_name_eff);
Bave=adapt_rda(Bave,NTIME,ndimrho,itotransposeback);
[B2ave,e]=rdaD3D_eff(shot,DIAG,'B2ave',exp_name_eff);
B2ave=adapt_rda(B2ave,NTIME,ndimrho,itotransposeback);
if strcmp(DIAG,'IDE')
FTRA.value=[];
else
[FTRA,e]=rdaD3D_eff(shot,DIAG,'FTRA',exp_name_eff);
FTRA=adapt_rda(FTRA,NTIME,ndimrho,itotransposeback);
end
else
Rinv.value=[]; R2inv.value=[]; Bave.value=[]; B2ave.value=[]; FTRA.value=[];
end
[LPFx,e]=rdaD3D_eff(shot,DIAG,'LPFx',exp_name_eff);
LPFx.value=LPFx.value(1:NTIME); LPFx.data=LPFx.value; LPFx.t=LPFx.t(1:NTIME);
[PFxx,e]=rdaD3D_eff(shot,DIAG,'PFxx',exp_name_eff);
PFxx=adapt_rda(PFxx,NTIME,max(LPFx.value)+1,itotransposeback);
[RPFx,e]=rdaD3D_eff(shot,DIAG,'RPFx',exp_name_eff);
RPFx=adapt_rda(RPFx,NTIME,max(LPFx.value)+1,itotransposeback);
[zPFx,e]=rdaD3D_eff(shot,DIAG,'zPFx',exp_name_eff);
zPFx=adapt_rda(zPFx,NTIME,max(LPFx.value)+1,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);
% 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
% set NaNs to zeroes
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
% get x values
gdat_data.psi(:,it)=psi_tree.value(it,Lpf1:-1:1)';
gdat_data.psi_axis(it)= gdat_data.psi(1,it);
gdat_data.psi_lcfs(it)= gdat_data.psi(end,it);
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(:,it) = sqrt(abs((gdat_data.psi(:,it)-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it))));
if strcmp(DIAG,'EQR'); gdat_data.rhopolnorm(end,it) = 1.;
% q value has only a few values and from center to edge, assume they are from central rhopol values on gdat_data.vol(:,it)=gdat_data.rhovolnorm(:,it).^2 .* eq.results.aeqdsk.volume(it);
% But they are every other point starting from 3rd gdat_data.Rmesh(:,it) = eq.results.geqdsk.r;
ijk=find(gdat_data.qvalue(:,it)~=0); gdat_data.Zmesh(:,it) = eq.results.geqdsk.z;
if length(ijk)>2 gdat_data.psi2D(:,:,it) = equil_all_t.gdata(it).psirz;
% now shots have non-zero axis values in eqr gdat_data.pressure(:,it) = equil_all_t.gdata(it).pres;
rhoeff=gdat_data.rhopolnorm(ijk,it); gdat_data.dpressuredpsi(:,it) = equil_all_t.gdata(it).pprime; % 2nd index are dV/dpsi
qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above gdat_data.ffprime(:,it) = equil_all_t.gdata(it).ffprim;
if ijk(1)>1 gdat_data.rwall(:,it) = equil_all_t.gdata(it).xlim;
rhoeff = [0.; rhoeff]; gdat_data.zwall(:,it) = equil_all_t.gdata(it).ylim;
qeff = [qeff(1) ;qeff]; if equil_all_t.adata(it).rxpt1 < -9.5
end gdat_data.R_Xpoints(:,it,1) = NaN;
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
gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)';
gdat_data.rhotornorm(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.phi(end,it)));
% get rhovol values
gdat_data.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)';
gdat_data.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
gdat_data.rhovolnorm(:,it) = sqrt(abs(gdat_data.vol(:,it) ./ gdat_data.vol(end,it)));
gdat_data.area(:,it)=Area.value(it,2*Lpf1-1:-2:1)';
gdat_data.dareadpsi(:,it)=Area.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
gdat_data.Rmesh(:,it) = Ri.value(it,1:M_Rmesh);
gdat_data.Zmesh(:,it) = Zj.value(it,1:N_Zmesh);
gdat_data.psi2D(1:M_Rmesh,1:N_Zmesh,it) = PFM_tree.value(1:M_Rmesh,1:N_Zmesh,it);
gdat_data.pressure(:,it)=Pres.value(it,2*Lpf1-1:-2:1)';
gdat_data.dpressuredpsi(:,it)=Pres.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
if ~isempty(Jpol.value)
gdat_data.jpol(:,it)=Jpol.value(it,2*Lpf1-1:-2:1)';
gdat_data.djpolpsi(:,it)=Jpol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
else
gdat_data.jpol = [];
gdat_data.djpolpsi = [];
end
gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)';
gdat_data.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1);
gdat_data.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1);
gdat_data.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1);
if ~isempty(Rinv.value)
gdat_data.rinv(:,it) = Rinv.value(it,Lpf1:-1:1)';
else
gdat_data.rinv = [];
end
if ~isempty(R2inv.value)
gdat_data.r2inv(:,it) = R2inv.value(it,Lpf1:-1:1)';
else else
gdat_data.r2inv = []; gdat_data.R_Xpoints(:,it,1) = equil_all_t.adata(it).rxpt1;
end end
if ~isempty(Bave.value) if equil_all_t.adata(it).rxpt1 < -9.5
gdat_data.bave(:,it) = Bave.value(it,Lpf1:-1:1)'; gdat_data.Z_Xpoints(:,it,1) = NaN;
else else
gdat_data.bave = []; gdat_data.Z_Xpoints(:,it,1) = equil_all_t.adata(it).zxpt1;
end end
if ~isempty(B2ave.value) if equil_all_t.adata(it).rxpt2 < -9.5
gdat_data.b2ave(:,it) = B2ave.value(it,Lpf1:-1:1)'; gdat_data.R_Xpoints(:,it,2) = NaN;
else else
gdat_data.b2ave = []; gdat_data.R_Xpoints(:,it,2) = equil_all_t.adata(it).rxpt2;
end end
if ~isempty(FTRA.value) if equil_all_t.adata(it).rxpt2 < -9.5
gdat_data.ftra(:,it) = FTRA.value(it,Lpf1:-1:1)'; gdat_data.Z_Xpoints(:,it,2) = NaN;
else else
gdat_data.ftra = []; gdat_data.Z_Xpoints(:,it,2) = equil_all_t.adata(it).zxpt2;
end end
% %
end end
gdat_data.x = gdat_data.rhopolnorm; gdat_data.x = gdat_data.rhopolnorm;
% %
[equil_Rcoil,e]=rdaD3D_eff(shot,DIAG,'Rcl',exp_name_eff); gdat_data.ip = [equil_all_t.gdata(:).cpasma];
gdat_data.Rcoils=equil_Rcoil.value;
[equil_Zcoil,e]=rdaD3D_eff(shot,DIAG,'Zcl',exp_name_eff);
gdat_data.Zcoils=equil_Zcoil.value;
%
if strcmp(DIAG,'IDE')
[IpiPSI,e]=rdaD3D_eff(shot,'IDG','Itor',exp_name_eff);
if length(IpiPSI.value)~=NTIME
disp('Problem with nb time points of IDG/Itor with respect to IDE NTIME, check with O. Sauter')
return
end
else
[IpiPSI,e]=rdaD3D_eff(shot,DIAG,'IpiPSI',exp_name_eff);
end
gdat_data.Ip = IpiPSI.value(1:NTIME);
% %
gdat_data.data = gdat_data.qvalue; % put q in data gdat_data.data = gdat_data.qvalue; % put q in data
gdat_data.units=[]; % not applicable gdat_data.units='q'; % not applicable
gdat_data.data_fullpath = [DIAG ' from expname: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)']; gdat_data.data_fullpath = ...
gdat_data.cocos = 17; % should check FPP ['q(:,time) in .data, all from [equil_all_t,neq,eq,ier]=read_mds_eq_func(shot,gdat_data.gdat_params.source,''d3d'');'];
gdat_data.cocos = 2;
gdat_data.dim{1} = gdat_data.x; gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t; gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits{1} = 'rhopolnorm'; gdat_data.dimunits{1} = 'rhopolnorm';
...@@ -1236,8 +1064,10 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') ...@@ -1236,8 +1064,10 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
% LOAD MULTI CHANNEL DATA ECS % LOAD MULTI CHANNEL DATA ECS
% powers, frequencies, etc % powers, frequencies, etc
params_eff = gdat_data.gdat_params; params_eff = gdat_data.gdat_params;
params_eff.data_request={'ECS','PECRH'}; params_eff.data_request={'rf' '\echpwrc'};
% pgyro tot in index=9 gyro_names={'leia','luke','scarecrow','tinman','chewbacca','nasa'};
power_names={'ecleifpwrc','eclukfpwrc','ecscafpwrc','ectinfpwrc','ecchefpwrc','ecnasfpwrc'};
% pgyro tot in index=length(gyro_names)+1
try try
gdat_data=gdat_d3d(shot,params_eff); gdat_data=gdat_d3d(shot,params_eff);
gdat_data.data_request = data_request_eff; gdat_data.data_request = data_request_eff;
...@@ -1248,129 +1078,75 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') ...@@ -1248,129 +1078,75 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
return return
end end
nb_timepoints = length(gdat_data.t); nb_timepoints = length(gdat_data.t);
pgyro = NaN*ones(nb_timepoints,9); pgyro = NaN*ones(nb_timepoints,length(gyro_names)+1);
pgyro(:,9) = reshape(gdat_data.data,nb_timepoints,1); pgyro(:,length(gyro_names)+1) = reshape(gdat_data.data,nb_timepoints,1);
gdat_data.data = pgyro; gdat_data.data = pgyro;
labels{9} = 'ECtot'; labels{length(gyro_names)+1} = 'ECtot';
for i=1:4 for i=1:length(gyro_names)
% "old" ECRH1 gyrotrons: gyro 1 to 4 in pgyro % "old" ECRH1 gyrotrons: gyro 1 to 4 in pgyro
params_eff.data_request={'ECS',['PG' num2str(i)]}; params_eff.data_request={'rf',['ech.' gyro_names{i} ':' power_names{i}]};
gdat_data_i=gdat_d3d(shot,params_eff); gdat_data_i=gdat_d3d(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else else
gdat_data.ec{i} = gdat_data_i; gdat_data.ec{i} = gdat_data_i;
gdat_data.data(:,i) = reshape(gdat_data_i.data,nb_timepoints,1); gdat_data.data(:,i) = reshape(gdat_data_i.data,nb_timepoints,1);
gdat_data.dim = [{gdat_data_i.t} {[1:9]}]; gdat_data.dim = [{gdat_data_i.t} {[1:length(gyro_names)+1]}];
if max(gdat_data_i.data) > 0. if max(gdat_data_i.data) > 0.
labels{i} = ['EC_' num2str(i)]; % labels{i} = [gyro_names{i} ':' power_names{i}];
labels{i} = [gyro_names{i}];
end end
end end
try try
a = sf2par('ECS',shot,'gyr_freq',['P_sy1_g' num2str(i)]); params_eff.data_request={'rf',['ech.' gyro_names{i} ':ec' gyro_names{i}(1:3) 'polang']};
a=gdat_d3d(shot,params_eff);
catch catch
% gyr_freq not present (old shots for example) % polang not present
a=[]; a=[];
end end
if isempty(a) if isempty(a)
else else
gdat_data.ec{i}.freq = a; gdat_data.ec{i}.polang_t = a.t;
gdat_data.freq_ec(i) = a.value; gdat_data.ec{i}.polang = a.data;
gdat_data.polang_ec{i} = a;
end end
try try
a = sf2par('ECS',shot,'GPolPos',['P_sy1_g' num2str(i)]); params_eff.data_request={'rf',['ech.' gyro_names{i} ':ec' gyro_names{i}(1:3) 'polcnt']};
a=gdat_d3d(shot,params_eff);
catch catch
% GPolPos not present
a=[]; a=[];
end end
if isempty(a) if isempty(a)
else else
gdat_data.ec{i}.polpos = a.value; gdat_data.ec{i}.polcnt = a.data;
gdat_data.polpos_ec(i) = a.value; gdat_data.polcnt_ec{i} = a;
end end
try try
a = sf2par('ECS',shot,'GTorPos',['P_sy1_g' num2str(i)]); params_eff.data_request={'rf',['ech.' gyro_names{i} ':ec' gyro_names{i}(1:3) 'torcnt']};
a=gdat_d3d(shot,params_eff);
catch catch
a=[]; a=[];
end end
if isempty(a) if isempty(a)
else else
gdat_data.ec{i}.torpos = a.value; gdat_data.ec{i}.torcnt = a.data;
gdat_data.torpos_ec(i) = a.value; gdat_data.torcnt_ec{i} = a;
end
% "new" ECRH2 gyrotrons: gyro 5 to 8 in pgyro
params_eff.data_request={'ECS',['PG' num2str(i) 'N']};
gdat_data_i=gdat_d3d(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i+4} = gdat_data_i;
gdat_data.data(:,i+4) = reshape(gdat_data_i.data,nb_timepoints,1);
gdat_data.dim = [{gdat_data_i.t} {[1:9]}];
if max(gdat_data_i.data) > 0.
labels{i+4} = ['EC_' num2str(i+4)];
end
end
try
a = sf2par('ECS',shot,'gyr_freq',['P_sy2_g' num2str(i)]);
catch
a=[];
end
if isempty(a)
else
gdat_data.ec{i+4}.freq = a;
gdat_data.freq_ec(i+4) = a.value;
end
try
a = sf2par('ECS',shot,'GPolPos',['P_sy2_g' num2str(i)]);
catch
a=[];
end
if isempty(a)
else
gdat_data.ec{i+4}.polpos = a.value;
gdat_data.polpos_ec(i+4) = a.value;
end
try
a = sf2par('ECS',shot,'GTorPos',['P_sy2_g' num2str(i)]);
catch
a=[];
end
if isempty(a)
else
gdat_data.ec{i+4}.torpos = a.value;
gdat_data.torpos_ec(i+4) = a.value;
end
params_eff.data_request={'ECN',['G' num2str(i) 'POL']};
gdat_data_i=gdat_d3d(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i+4}.gpol_ec = gdat_data_i;
end
params_eff.data_request={'ECN',['G' num2str(i) 'TOR']};
gdat_data_i=gdat_d3d(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i+4}.gtor_ec = gdat_data_i;
end
params_eff.data_request={'ECN',['G' num2str(i) 'PO4']};
gdat_data_i=gdat_d3d(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i+4}.gpo4_ec = gdat_data_i;
end
params_eff.data_request={'ECN',['G' num2str(i) 'PO8']};
gdat_data_i=gdat_d3d(shot,params_eff);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i+4}.gpo8_ec = gdat_data_i;
end end
end end
gdat_data.gyro_names = gyro_names;
gdat_data.power_names = power_names;
% add ech on time_interval from total power>3% of max;
ij=find(gdat_data.data(:,end)>0.03.*max(gdat_data.data(:,end)));
if ~isempty(ij)
gdat_data.ec_t_on = [gdat_data.t(ij(1)) gdat_data.t(ij(end))];
else
gdat_data.ec_t_on = [0 0];
end
if ~isempty(gdat_data.dim) if ~isempty(gdat_data.dim)
gdat_data.t = gdat_data.dim{1}; gdat_data.t = gdat_data.dim{1};
gdat_data.x = gdat_data.dim{2}; gdat_data.x = gdat_data.dim{2};
gdat_data.dimunits=[{'time [s]'} {'ECRH1(1:4) ECRH2(1:4) ECtot'}]; gdat_data.dimunits=[{'time [ms]'} gyro_names{:} {'ECtot'}];
gdat_data.units='W'; gdat_data.units='W';
gdat_data.freq_ech_units = 'GHz'; gdat_data.data_fullpath=['rf::ech.gyro_names:ec..pwrc'];
gdat_data.data_fullpath=['ECS/' 'PGi and PGiN, etc'];
icount=0; icount=0;
for i=1:length(labels) for i=1:length(labels)
if ~isempty(labels{i}) if ~isempty(labels{i})
...@@ -1381,7 +1157,6 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') ...@@ -1381,7 +1157,6 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
else else
gdat_data.freq_ech_units =[]'; gdat_data.freq_ech_units =[]';
gdat_data.ec = []; gdat_data.ec = [];
gdat_data.freq_ec = [];
gdat_data.polpos_ec = []; gdat_data.polpos_ec = [];
gdat_data.torpos_ec = []; gdat_data.torpos_ec = [];
end end
...@@ -1756,67 +1531,22 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') ...@@ -1756,67 +1531,22 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'sxr'} case {'sxr'}
% sxr from B by default or else if 'source','else' is provided % sxr from sx90rm1s by default or else if 'source' is provided
if ~isfield(gdat_data.gdat_params,'freq')|| isempty(gdat_data.gdat_params.freq) if ~isfield(gdat_data.gdat_params,'freq')|| isempty(gdat_data.gdat_params.freq)
gdat_data.gdat_params.freq = 'slow'; gdat_data.gdat_params.freq = 1;
end end
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = 'G'; gdat_data.gdat_params.source = 'sx90rm1s';
elseif length(gdat_data.gdat_params.source)>1 || ...
upper(gdat_data.gdat_params.source)<'F' || upper(gdat_data.gdat_params.source)>'M'
if gdat_data.gdat_params.nverbose>=1
warning(['source = ' gdat_data.gdat_params.source ' not expected with data_request= ' data_request_eff])
end
return
end end
dim1_len = 'from_chord_nb'; % thus up to max(chords_ok_nb)
if ~isfield(gdat_data.gdat_params,'camera') || isempty(gdat_data.gdat_params.camera) if ~isfield(gdat_data.gdat_params,'camera') || isempty(gdat_data.gdat_params.camera)
gdat_data.gdat_params.camera = []; gdat_data.gdat_params.camera = [1:28];
elseif ischar(gdat_data.gdat_params.camera) && strcmp(gdat_data.gdat_params.camera,'central')
gdat_data.gdat_params.source = 'J';
gdat_data.gdat_params.camera = 49;
elseif isnumeric(gdat_data.gdat_params.camera)
% ok keep the array, but first dim to contain just the related chords
dim1_len='nb_of_chords';
else
if gdat_data.gdat_params.nverbose>=1
warning(['camera = ' gdat_data.gdat_params.camera ' not expected with data_request= ' data_request_eff])
end
return
end end
if length(gdat_data.gdat_params.camera)==1; dim1_len='nb_of_chords'; end
gdat_data.gdat_params.source = upper(gdat_data.gdat_params.source); gdat_data.gdat_params.source = upper(gdat_data.gdat_params.source);
% %
if ~isfield(gdat_data.gdat_params,'time_interval') if ~isfield(gdat_data.gdat_params,'time_interval')
gdat_data.gdat_params.time_interval = []; gdat_data.gdat_params.time_interval = [];
end end
[aa,bb]=unix(['ssh ' '-o "StrictHostKeyChecking no" sxd3d20.d3d.ipp.mpg.de WhichSX ' num2str(shot) ' ' ... exp_name_eff = 'D3D';
upper(gdat_data.gdat_params.source)]);
chords_ok_diags=regexpi(bb,'(?<chord>[A-Z]_[0-9]+)\saddress:\s[0-9]+\sDiag:\s(?<diag>[A-Z]+)','names');
if isempty(chords_ok_diags)
chords_ok_nb = [];
else
for i=1:length(chords_ok_diags)
chords_ok_nb(i) = str2num(chords_ok_diags(i).chord(3:end));
end
end
if isempty(gdat_data.gdat_params.camera);
gdat_data.gdat_params.camera = chords_ok_nb;
else
for i=1:length(gdat_data.gdat_params.camera)
ij=find(chords_ok_nb==gdat_data.gdat_params.camera(i));
if ~isempty(ij)
chords_ok_diags_tmp(i) = chords_ok_diags(ij);
else
% chord not in use
chords_ok_diags_tmp(i).chord = [];
chords_ok_diags_tmp(i).diag = [];
end
end
chords_ok_diags = chords_ok_diags_tmp;
chords_ok_nb = gdat_data.gdat_params.camera;
end
exp_name_eff = 'D3DD';
icount = 0; icount = 0;
nnth = 1; nnth = 1;
if isnumeric(gdat_data.gdat_params.freq) && gdat_data.gdat_params.freq>1; if isnumeric(gdat_data.gdat_params.freq) && gdat_data.gdat_params.freq>1;
...@@ -1824,56 +1554,34 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') ...@@ -1824,56 +1554,34 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
gdat_data.gdat_params.freq = nnth; gdat_data.gdat_params.freq = nnth;
end end
for i=1:length(gdat_data.gdat_params.camera) for i=1:length(gdat_data.gdat_params.camera)
if ischar(gdat_data.gdat_params.freq) && strcmp(gdat_data.gdat_params.freq,'slow'); chords_ok_diags(i).diag = 'SSX'; end tree = 'spectroscopy';
if ~isempty(chords_ok_diags(i).diag) && ~isempty(chords_ok_diags(i).chord) ichord = gdat_data.gdat_params.camera(i);
[a,e]=rdaD3D_eff(shot,chords_ok_diags(i).diag,chords_ok_diags(i).chord,exp_name_eff,gdat_data.gdat_params.time_interval); diagname = ['sxr:' gdat_data.gdat_params.source ':' gdat_data.gdat_params.source num2str(ichord,'%.2d')];
else a = gdat_d3d(shot,{tree,diagname});
a.data = [];
a.t = [];
end
if ~isempty(a.data) if ~isempty(a.data)
icount = icount + 1; icount = icount + 1;
if icount == 1 if icount == 1
% first time has data % first time has data
% assume all chords have same time base even if from different shotfile
% time missing one point
if length(a.value) == length(a.t)+1
a.t=linspace(a.range(1),a.range(2),length(a.value));
a_time.index(2) = length(a.value);
end
gdat_data.t = a.t(1:nnth:end); gdat_data.t = a.t(1:nnth:end);
gdat_data.units = a.units; gdat_data.units = a.units;
if strcmp(dim1_len,'from_chord_nb') gdat_data.data = NaN*ones(max(gdat_data.gdat_params.camera),length(gdat_data.t));
gdat_data.data = NaN*ones(max(chords_ok_nb),length(gdat_data.t)); gdat_data.x = [1:max(gdat_data.gdat_params.camera)]; % simpler to have index corresponding to chord number, except for central
gdat_data.dim{1} = [1:max(chords_ok_nb)]; % simpler to have index corresponding to chord number, except for central gdat_data.dim{1} = gdat_data.x;
else
gdat_data.data = NaN*ones(length(chords_ok_nb),length(gdat_data.t));
gdat_data.dim{1} = chords_ok_nb;
end
gdat_data.dim{2} = gdat_data.t; gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits = [{'chord nb'}; {'s'}]; gdat_data.dimunits = [{'chord nb'}; {'s'}];
gdat_data.data_fullpath = ['sxr from source=''' gdat_data.gdat_params.source ' uses shotfiles SX' ... gdat_data.data_fullpath = ['sxr from source=''' gdat_data.gdat_params.source ''''];
setdiff(unique(strvcat(chords_ok_diags(:).diag)),'SX')'];
gdat_data.label = ['SXR/' upper(gdat_data.gdat_params.source)]; gdat_data.label = ['SXR/' upper(gdat_data.gdat_params.source)];
end end
try try
if strcmp(dim1_len,'from_chord_nb') gdat_data.data(ichord,:) = a.data(1:nnth:end);
gdat_data.data(chords_ok_nb(i),:) = a.data(1:nnth:end);
else
gdat_data.data(i,:) = a.data(1:nnth:end);
end
catch catch
if (gdat_params.nverbose>=1); disp(['problem with ' chords_ok_diags(i).diag ' ; ' chords_ok_diags(i).chord]); end if (gdat_params.nverbose>=1); disp(['problem with ' gdat_data.gdat_params.source '...' num2str(ichord)]); end
end end
else else
% add fields not yet defined in case all cases have empty data % add fields not yet defined in case all cases have empty data
end end
end end
gdat_data.chords = chords_ok_nb; gdat_data.chords = gdat_data.gdat_params.camera;
if length(gdat_data.dim)>=1
gdat_data.x = gdat_data.dim{1};
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'transp'} case {'transp'}
...@@ -1936,36 +1644,3 @@ return ...@@ -1936,36 +1644,3 @@ return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% functions for portions called several times % functions for portions called several times
function [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data);
% get basic params to be able to read results in EQ-like shotfiles
% M_Rmesh,N_Zmesh only needed for equil when 2D quantities are required
%
shot=gdat_data.shot;
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 = rdaD3D_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_data.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_data.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
[equil_time,e]=rdaD3D_eff(shot,DIAG,'time',exp_name_eff);
gdat_data.t = equil_time.value(1:NTIME);
crpptbx/D3D/get_eqdsk_d3d.m 0 → 100644
+ 128
− 0
View file @ cde58f23
function [eqdskd3d, equil_all_t, equil_t_index]=get_eqdsk_d3d(shot,time,zshift,varargin);
%
% called by:
%
% eqdsk = gdat(shot/equil,'eqdsk','time',time[,'source','efit02',...]);
%
% [eqdskd3d, equil_all_t, equil_t_index]=geteqdskd3d(shot,time,zshift,varargin);
%
% if shot is a structure assumes obtained from gdat(shot,'equil',...);
%
% varargin{1:2}: 'source','efit03' (default) or 'source','efit01' (case insensitive)
%
% [eqdskd3d, equil_all_t, equil_t_index]=geteqdskd3d(shot,time);
% [eqdskd3d, equil_all_t, equil_t_index]=geteqdskd3d(shot,time,[],'source','efit02');
%
% you can then do:
% write_eqdsk('fname',eqdskd3d,2); % EFIT is COCOS=2 by default (except abs(q) is set, so now q has a sign)
% write_eqdsk('fname',eqdskd3d,[2 11]); % if you want an ITER version with COCOS=11
%
if ~exist('shot') || isempty(shot);
disp('requires a shot or equil structure from gdat(shot,''equil'',...) in geteqdskd3d')
return
end
if ~exist('time');
time_eff = 2.0*1e3;
else
time_eff = time;
end
if ~exist('zshift') || isempty(zshift) || ~isnumeric(zshift)
zshift_eff = 0.; % no shift
else
zshift_eff = zshift;
end
if nargin >= 5 && ~isempty(varargin{1})
if ~isempty(varargin{2})
equil_source = varargin{2};
else
disp(['Warning source for equil in get_eqdsk_d3d not defined']);
return;
end
else
equil_source = 'efit03';
end
if isnumeric(shot)
% check read_mds_eq_func is available
equil_in = gdat(shot,'equil','equil',equil_source);
else
equil_in = shot;
shot = equil_in.shot;
end
equil_all_t = equil_in.equil_all_t;
equil = equil_all_t;
[zz it]=min(abs(equil.time-time_eff));
equil_t_index = it;
eqdsk.cocos=2;
eqdsk.nr = equil.gdata(it).nw;
eqdsk.nz = equil.gdata(it).nh;
eqdsk.rboxlen = equil.gdata(it).xdim ;
eqdsk.rboxleft = equil.gdata(it).rgrid1;
eqdsk.zboxlen = equil.gdata(it).zdim ;
eqdsk.zmid = equil.gdata(it).zmid ;
eqdsk.rmesh = linspace(eqdsk.rboxleft,eqdsk.rboxleft+eqdsk.rboxlen,eqdsk.nr)';
eqdsk.zmesh = linspace(eqdsk.zmid-eqdsk.zboxlen/2,eqdsk.zmid+eqdsk.zboxlen/2,eqdsk.nz)' - zshift_eff;
eqdsk.zshift = zshift_eff;
eqdsk.psi = equil.gdata(it).psirz;
eqdsk.psirz = reshape(eqdsk.psi,prod(size(eqdsk.psi)),1);
eqdsk.psiaxis = equil.gdata(it).ssimag;
eqdsk.psiedge = equil.gdata(it).ssibry;
% need psimesh ot be equidistant and same nb of points as nr
eqdsk.psimesh=linspace(0,1,eqdsk.nr)';
eqdsk.rhopsi = sqrt(eqdsk.psimesh);
% psimesh assumed from psi_axis on axis to have correct d/dpsi
eqdsk.psimesh = (eqdsk.psiedge-eqdsk.psiaxis) .* eqdsk.psimesh + eqdsk.psiaxis;
eqdsk.p = equil.gdata(it).pres;
eqdsk.pprime = - equil.gdata(it).pprime; % pprime has wrong sign
eqdsk.FFprime = - equil.gdata(it).ffprim; % ffprime has wrong sign
eqdsk.F = equil.gdata(it).fpol;
eqdsk.q = equil.gdata(it).qpsi.*sign(equil.gdata(it).bzero).*sign(equil.gdata(it).cpasma); % EFIT uses always abs(q) so recover sign
eqdsk.b0 = equil.gdata(it).bzero;
eqdsk.r0 = equil.gdata(it).rzero;
eqdsk.ip = equil.gdata(it).cpasma;
eqdsk.raxis = equil.gdata(it).rmaxis;
eqdsk.zaxis = equil.gdata(it).zmaxis - eqdsk.zshift;
eqdsk.nbbound = equil.gdata(it).nbbbs;
eqdsk.rplas = equil.gdata(it).rbbbs(1:equil.gdata(it).nbbbs);
eqdsk.zplas = equil.gdata(it).zbbbs(1:equil.gdata(it).nbbbs);
eqdsk.nblim = equil.gdata(it).limitr;
eqdsk.rlim = equil.gdata(it).xlim(1:equil.gdata(it).limitr);
eqdsk.zlim = equil.gdata(it).ylim(1:equil.gdata(it).limitr);
eqdsk.stitle=['#' num2str(shot) ' t=' num2str(time_eff) ' from ' equil_in.gdat_params.exp_name '/' equil_in.gdat_params.equil];
eqdsk.ind1=3;
psisign = sign(eqdsk.psimesh(end)-eqdsk.psimesh(1));
[dum1,dum2,dum3,F2_05]=interpos(psisign.*eqdsk.psimesh,eqdsk.FFprime,-0.1);
fedge=eqdsk.r0.*eqdsk.b0;
F2 = psisign.*2.*F2_05 + fedge.^2;
if abs(eqdsk.F(end)-fedge)./abs(fedge)>1e-6
disp('problem with cocos?')
disp(['eqdsk.F(end) = ' num2str(eqdsk.F)])
disp(['r0*b0 = ' num2str(eqdsk.r0) '*' num2str(eqdsk.b0) ' = ' num2str(fedge)]);
end
% get plasma boundary
figure
contour(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',100)
hold all
plot(eqdsk.rplas,eqdsk.zplas)
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
plot(eqdsk.rlim,eqdsk.zlim,'k-')
title(eqdsk.stitle)
eqdskd3d = eqdsk;
crpptbx/D3D/get_signal_d3d.m
+ 18
− 17
View file @ cde58f23
...@@ -3,26 +3,27 @@ function [gdat_data,error_status] = get_signal_d3d(signal,gdat_data); ...@@ -3,26 +3,27 @@ function [gdat_data,error_status] = get_signal_d3d(signal,gdat_data);
gdat_data.data = mdsvalue(signal); gdat_data.data = mdsvalue(signal);
error_status = 1; error_status = 1;
gdat_data.dim = [];
gdat_data.dimunits = [];
gdat_data.t = [];
gdat_data.x = [];
nbdims=numel(setdiff(size(gdat_data.data),1)); nbdims=numel(setdiff(size(gdat_data.data),1));
if ~isempty(gdat_data.data) && ~ischar(gdat_data.data) if ~isempty(gdat_data.data) && ~ischar(gdat_data.data)
gdat_data.units = mdsvalue(['units_of(' signal ')']); gdat_data.units = mdsvalue(['units_of(' signal ')']);
for i=1:nbdims if prod(size(gdat_data.data)) > 1
gdat_data.dim{i} = mdsvalue(['dim_of(' signal ',' num2str(i-1) ')']); for i=1:nbdims
% gdat_data.dimunits{i} = mdsvalue(['dimunits_of(' signal ',' num2str(i-1) ')']); gdat_data.dim{i} = mdsvalue(['dim_of(' signal ',' num2str(i-1) ')']);
gdat_data.dimunits{i} = 'to get'; % gdat_data.dimunits{i} = mdsvalue(['dimunits_of(' signal ',' num2str(i-1) ')']);
gdat_data.dimunits{i} = 'to get';
end
if nbdims==1
gdat_data.t = gdat_data.dim{1};
gdat_data.x = [];
else
gdat_data.t = gdat_data.dim{1};
gdat_data.x = gdat_data.dim{2};
end
error_status = 0;
end end
if nbdims==1
gdat_data.t = gdat_data.dim{1};
gdat_data.x = [];
else
gdat_data.t = gdat_data.dim{1};
gdat_data.x = gdat_data.dim{2};
end
error_status = 0;
else
gdat_data.dim = [];
gdat_data.dimunits = [];
gdat_data.t = [];
gdat_data.x = [];
end end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment