diff --git a/crpptbx/D3D/d3d_requests_mapping.m b/crpptbx/D3D/d3d_requests_mapping.m
index 71e14c8630f7797976dbf883785b046ede608b3b..4ddb0741403e494f2ca3efdc69bd71aad5916abc 100644
--- a/crpptbx/D3D/d3d_requests_mapping.m
+++ b/crpptbx/D3D/d3d_requests_mapping.m
@@ -52,7 +52,7 @@ switch lower(data_request)
mapping.timedim = 1;
mapping.label = 'B_0';
mapping.method = 'signal';
- mapping.expression = [{'FPC'},{'BTF'}];
+ mapping.expression = [{'EFIT01'},{'\bcentr'}];
case 'beta'
mapping.timedim = 1;
mapping.label = '\beta';
@@ -106,23 +106,23 @@ switch lower(data_request)
mapping.label = 'delta\_top';
mapping.timedim = 1;
mapping.method = 'signal';
- mapping.expression = [{'FPG'},{'delRoben'}];
+ mapping.expression = [{'OT01'},{'\tritop'}];
case 'delta_bottom'
mapping.label = 'delta\_bottom';
mapping.timedim = 1;
mapping.method = 'signal';
- mapping.expression = [{'FPG'},{'delRuntn'}];
- case {'ece', 'eced', 'ece_rho', 'eced_rho'}
- mapping.timedim = 2;
+ mapping.expression = [{'OT01'},{'\tribot'}];
+ case {'ece', 'ece_rho'}
+ mapping.timedim = 1;
mapping.method = 'switchcase';
mapping.expression = '';
case 'eqdsk'
mapping.timedim = 2;
- mapping.method = 'switchcase'; % could use function make_eqdsk directly?
+ mapping.method = 'switchcase';
mapping.expression = '';
case 'equil'
mapping.gdat_timedim = 2;
- mapping.method = 'switchcase'; % could use function make_eqdsk directly?
+ mapping.method = 'switchcase';
mapping.expression = '';
case 'halpha'
mapping.timedim = 1;
@@ -148,12 +148,12 @@ switch lower(data_request)
mapping.timedim = 1;
mapping.label = 'Plasma current';
mapping.method = 'signal';
- mapping.expression = [{'MAG'},{'Ipa'}];
+ mapping.expression = [{'efit01'},{'\cpasma'}];
case 'kappa'
mapping.timedim = 1;
mapping.label = '\kappa';
mapping.method = 'signal';
- mapping.expression = [{'FPG'},{'k'}];
+ mapping.expression = [{'efit01'},{'\kappa'}];
case 'kappa_top'
mapping.timedim = 1;
mapping.label = '\kappa^{top}';
@@ -168,7 +168,7 @@ switch lower(data_request)
mapping.timedim = 1;
mapping.label = 'l_i';
mapping.method = 'signal';
- mapping.expression = [{'FPG'},{'li'}];
+ mapping.expression = [{'efit01'},{'\li3'}];
case 'mhd'
mapping.timedim = 1;
mapping.label = 'Odd and Even n';
@@ -188,6 +188,11 @@ switch lower(data_request)
'gdat_tmp.n_even.amp_t=gdat_tmp2.t;' ...
'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;'];
+ case 'nbi'
+ mapping.timedim = 1;
+ mapping.label = 'Pnbi';
+ mapping.method = 'signal';
+ mapping.expression = [{'nb'},{'pinj'}];
case 'ne'
mapping.timedim = 2;
mapping.method = 'switchcase';
@@ -198,14 +203,9 @@ switch lower(data_request)
mapping.expression = [{'DCN'},{'H-1'}];
case 'nel'
mapping.timedim = 1;
- mapping.label = 'line-averaged el. density';
- mapping.expression = [{'FPG'},{'lenH-1'}];
- mapping.method = 'expression';
- 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);'];
+ mapping.label = 'NEBAR\_R0';
+ mapping.method = 'signal';
+ mapping.expression = [{'efit01'},{'\NEBAR_R0'}];
case 'ne_rho'
mapping.timedim = 2;
mapping.label = 'ne';
@@ -254,6 +254,11 @@ switch lower(data_request)
mapping.gdat_timedim = 2;
mapping.label = 'q';
mapping.method = 'switchcase';
+ case 'r0'
+ mapping.timedim = 1;
+ mapping.label = 'R_0';
+ mapping.method = 'signal';
+ mapping.expression = [{'EFIT01'},{'\rzero'}];
case 'rgeom'
mapping.label = 'Rgeom';
mapping.timedim = 1;
diff --git a/crpptbx/D3D/gdat_d3d.m b/crpptbx/D3D/gdat_d3d.m
index 1b640907a3e531463dadb0f070c96d06d40da340..d20e5db508939d8ce75eabeec0d2b30b4057bdb1 100644
--- a/crpptbx/D3D/gdat_d3d.m
+++ b/crpptbx/D3D/gdat_d3d.m
@@ -198,9 +198,7 @@ if (nargin>=ivarargin_first_char)
end
data_request_eff = gdat_params.data_request; % in case was defined in pairs
% default equil:
-if ~isfield(gdat_params,'equil'); gdat_params.equil = 'EQI'; 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 ~isfield(gdat_params,'equil'); gdat_params.equil = 'EFIT01'; end
% if it is a request_keyword can obtain description:
if ischar(data_request_eff) || length(data_request_eff)==1
@@ -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
return
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;
% end of method "signal"
@@ -533,7 +533,7 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- case {'ece', 'eced', 'ece_rho', 'eced_rho'}
+ case {'ece', 'ece_rho'}
nth_points = 13;
if isfield(gdat_data.gdat_params,'nth_points') && ~isempty(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')
if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels)
channels = gdat_data.gdat_params.channels;
end
+ diag_name = 'ece';
if nth_points>=10
match_rz_to_time = 1;
else
@@ -561,126 +562,66 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
gdat_data.gdat_params.time_interval = time_interval;
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
- channels = [1:inb_chord];
+ channels = [1:40];
end
+ chanelles=sort(channels);
gdat_data.dim{1} = channels;
gdat_data.gdat_params.channels = channels;
- if nth_points>1
- gdat_data.data = a.data(channels,1:nth_points:end);
- gdat_data.dim{2} = a.t(1:nth_points:end);
- else
- gdat_data.data = a.data(channels,:);
- gdat_data.dim{2} = a.t;
+ for i=channels
+ a = gdat_d3d(shot,{'ece',['\tece' num2str(i,'%.2d')]});
+ if nth_points>1
+ gdat_data.data(i,:) = a.data(1:nth_points:end);
+ gdat_data.dim{2} = a.t(1:nth_points:end);
+ else
+ gdat_data.data(i,:) = a.data;
+ gdat_data.dim{2} = a.t;
+ end
end
gdat_data.x = gdat_data.dim{1};
gdat_data.t = gdat_data.dim{2};
gdat_data.dimunits=[{'channels'} ; {'time [s]'}];
if any(strcmp(fieldnames(a),'units')); gdat_data.units=a.units; end
- try
- [aR,e]=rdaD3D_eff(shot,diag_name,'R-A',exp_name_eff,time_interval);
- catch
- end
- try
- [aZ,e]=rdaD3D_eff(shot,diag_name,'z-A',exp_name_eff,time_interval);
- catch
- disp(['problem with getting z-A in ' diag_name])
- end
- if match_rz_to_time
- % interpolate R structure on ece data time array, to ease plot vs R
- for i=1:length(channels)
- 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
+% $$$ if match_rz_to_time
+% $$$ % interpolate R structure on ece data time array, to ease plot vs R
+% $$$ for i=1:length(channels)
+% $$$ 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);
gdat_data.data(ij)=NaN;
- gdat_data.rz.r=radius.data;
- gdat_data.rz.z=zheight.data;
- 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.rz.r=radius.data;
+% $$$ gdat_data.rz.z=zheight.data;
+% $$$ gdat_data.rz.t = radiuszheight.t;
+ 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')
- params_equil = gdat_data.gdat_params;
- params_equil.data_request = 'equil';
- [equil,params_equil,error_status] = gdat_d3d(shot,params_equil);
- if error_status>0
- if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); 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
+% $$$ 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'}
%
- time_eqdsks=1.; % default time
+ time_eqdsks=2000.; % default time
if isfield(gdat_data.gdat_params,'time') && ~isempty(gdat_data.gdat_params.time)
time_eqdsks = gdat_data.gdat_params.time;
else
gdat_data.gdat_params.time = time_eqdsks;
disp(['"time" is expected as an option, choose default time = ' num2str(time_eqdsks)]);
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.t = time_eqdsks;
zshift = 0.;
@@ -702,7 +643,7 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
for itime=1:length(time_eqdsks)
time_eff = time_eqdsks(itime);
% 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')]);
cocos_out = equil.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')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'equil'}
- % get equil params and time array in gdat_data.t
- [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data);
- % 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
+ if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || ~ischar(gdat_data.gdat_params.source)
+ gdat_data.gdat_params.source = 'efit03';
else
- itotransposeback = 0;
end
- qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback);
- ijnan=find(isnan(qpsi.value));
- qpsi.value(ijnan)=0;
- [psi_tree,e]=rdaD3D_eff(shot,DIAG,'PFL',exp_name_eff);
- psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
- [phi_tree,e]=rdaD3D_eff(shot,DIAG,'TFLx',exp_name_eff);
- phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
- [Vol,e]=rdaD3D_eff(shot,DIAG,'Vol',exp_name_eff);
- Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
- [Area,e]=rdaD3D_eff(shot,DIAG,'Area',exp_name_eff);
- Area=adapt_rda(Area,NTIME,2*ndimrho,itotransposeback);
- [Ri,e]=rdaD3D_eff(shot,DIAG,'Ri',exp_name_eff);
- Ri=adapt_rda(Ri,NTIME,M_Rmesh,itotransposeback);
- [Zj,e]=rdaD3D_eff(shot,DIAG,'Zj',exp_name_eff);
- Zj=adapt_rda(Zj,NTIME,N_Zmesh,itotransposeback);
- [PFM_tree,e]=rdaD3D_eff(shot,DIAG,'PFM','-raw','-exp',exp_name_eff); % -raw necessary for IDE
- PFM_tree=adaptPFM_rda(PFM_tree,M_Rmesh,N_Zmesh,NTIME);
- [Pres,e]=rdaD3D_eff(shot,DIAG,'Pres',exp_name_eff);
- 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);
+ if exist('read_mds_eq_func')<=0
+ if gdat_params.nverbose>=3; disp(['addpath(''/m/GAtools/matlab/efit'',''end'') since needs function read_mds_eq_func']); end
+ addpath('/m/GAtools/matlab/efit','end');
+ end
+ [equil_all_t,neq,eq,ier]=read_mds_eq_func(shot,gdat_data.gdat_params.source,'d3d');
+ gdat_data.rhovolnorm = eq.results.geqdsk.rhovn;
+ gdat_data.equil_all_t = equil_all_t;
+ gdat_data.eq_raw = eq;
+ % extract main parameters, can add more later
+ gdat_data.phi = [];
+ gdat_data.rhotornorm = [];
+ gdat_data.t = equil_all_t.time;
+ for it=1:length(equil_all_t.time)
+ gdat_data.qvalue(:,it) = equil_all_t.gdata(it).qpsi;
+ gdat_data.psi(:,it) = (equil_all_t.gdata(it).ssibry - equil_all_t.gdata(it).ssimag) ...
+ .* linspace(0,1,equil_all_t.gdata(it).nw)' + equil_all_t.gdata(it).ssimag;
+ 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))));
- 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
- 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)';
+ gdat_data.rhopolnorm(end,it) = 1.;
+ 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;
+ gdat_data.psi2D(:,:,it) = equil_all_t.gdata(it).psirz;
+ gdat_data.pressure(:,it) = equil_all_t.gdata(it).pres;
+ gdat_data.dpressuredpsi(:,it) = equil_all_t.gdata(it).pprime; % 2nd index are dV/dpsi
+ gdat_data.ffprime(:,it) = equil_all_t.gdata(it).ffprim;
+ gdat_data.rwall(:,it) = equil_all_t.gdata(it).xlim;
+ gdat_data.zwall(:,it) = equil_all_t.gdata(it).ylim;
+ if equil_all_t.adata(it).rxpt1 < -9.5
+ gdat_data.R_Xpoints(:,it,1) = NaN;
else
- gdat_data.r2inv = [];
+ gdat_data.R_Xpoints(:,it,1) = equil_all_t.adata(it).rxpt1;
end
- if ~isempty(Bave.value)
- gdat_data.bave(:,it) = Bave.value(it,Lpf1:-1:1)';
+ if equil_all_t.adata(it).rxpt1 < -9.5
+ gdat_data.Z_Xpoints(:,it,1) = NaN;
else
- gdat_data.bave = [];
+ gdat_data.Z_Xpoints(:,it,1) = equil_all_t.adata(it).zxpt1;
end
- if ~isempty(B2ave.value)
- gdat_data.b2ave(:,it) = B2ave.value(it,Lpf1:-1:1)';
+ if equil_all_t.adata(it).rxpt2 < -9.5
+ gdat_data.R_Xpoints(:,it,2) = NaN;
else
- gdat_data.b2ave = [];
+ gdat_data.R_Xpoints(:,it,2) = equil_all_t.adata(it).rxpt2;
end
- if ~isempty(FTRA.value)
- gdat_data.ftra(:,it) = FTRA.value(it,Lpf1:-1:1)';
+ if equil_all_t.adata(it).rxpt2 < -9.5
+ gdat_data.Z_Xpoints(:,it,2) = NaN;
else
- gdat_data.ftra = [];
+ gdat_data.Z_Xpoints(:,it,2) = equil_all_t.adata(it).zxpt2;
end
%
end
gdat_data.x = gdat_data.rhopolnorm;
%
- [equil_Rcoil,e]=rdaD3D_eff(shot,DIAG,'Rcl',exp_name_eff);
- 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.ip = [equil_all_t.gdata(:).cpasma];
%
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.units='q'; % not applicable
+ gdat_data.data_fullpath = ...
+ ['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{2} = gdat_data.t;
gdat_data.dimunits{1} = 'rhopolnorm';
@@ -1236,8 +1064,10 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
% LOAD MULTI CHANNEL DATA ECS
% powers, frequencies, etc
params_eff = gdat_data.gdat_params;
- params_eff.data_request={'ECS','PECRH'};
- % pgyro tot in index=9
+ params_eff.data_request={'rf' '\echpwrc'};
+ 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
gdat_data=gdat_d3d(shot,params_eff);
gdat_data.data_request = data_request_eff;
@@ -1248,129 +1078,75 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
return
end
nb_timepoints = length(gdat_data.t);
- pgyro = NaN*ones(nb_timepoints,9);
- pgyro(:,9) = reshape(gdat_data.data,nb_timepoints,1);
+ pgyro = NaN*ones(nb_timepoints,length(gyro_names)+1);
+ pgyro(:,length(gyro_names)+1) = reshape(gdat_data.data,nb_timepoints,1);
gdat_data.data = pgyro;
- labels{9} = 'ECtot';
- for i=1:4
+ labels{length(gyro_names)+1} = 'ECtot';
+ for i=1:length(gyro_names)
% "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);
if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim)
else
gdat_data.ec{i} = gdat_data_i;
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.
- labels{i} = ['EC_' num2str(i)];
+ % labels{i} = [gyro_names{i} ':' power_names{i}];
+ labels{i} = [gyro_names{i}];
end
end
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
- % gyr_freq not present (old shots for example)
+ % polang not present
a=[];
end
if isempty(a)
else
- gdat_data.ec{i}.freq = a;
- gdat_data.freq_ec(i) = a.value;
+ gdat_data.ec{i}.polang_t = a.t;
+ gdat_data.ec{i}.polang = a.data;
+ gdat_data.polang_ec{i} = a;
end
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
- % GPolPos not present
a=[];
end
if isempty(a)
else
- gdat_data.ec{i}.polpos = a.value;
- gdat_data.polpos_ec(i) = a.value;
+ gdat_data.ec{i}.polcnt = a.data;
+ gdat_data.polcnt_ec{i} = a;
end
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
a=[];
end
if isempty(a)
else
- gdat_data.ec{i}.torpos = a.value;
- gdat_data.torpos_ec(i) = a.value;
- 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;
+ gdat_data.ec{i}.torcnt = a.data;
+ gdat_data.torcnt_ec{i} = a;
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)
gdat_data.t = gdat_data.dim{1};
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.freq_ech_units = 'GHz';
- gdat_data.data_fullpath=['ECS/' 'PGi and PGiN, etc'];
+ gdat_data.data_fullpath=['rf::ech.gyro_names:ec..pwrc'];
icount=0;
for i=1:length(labels)
if ~isempty(labels{i})
@@ -1381,7 +1157,6 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
else
gdat_data.freq_ech_units =[]';
gdat_data.ec = [];
- gdat_data.freq_ec = [];
gdat_data.polpos_ec = [];
gdat_data.torpos_ec = [];
end
@@ -1756,67 +1531,22 @@ elseif strcmp(mapping_for_d3d.method,'switchcase')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)
- gdat_data.gdat_params.freq = 'slow';
+ gdat_data.gdat_params.freq = 1;
end
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
- gdat_data.gdat_params.source = 'G';
- 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
+ gdat_data.gdat_params.source = 'sx90rm1s';
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)
- gdat_data.gdat_params.camera = [];
- 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
+ gdat_data.gdat_params.camera = [1:28];
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);
%
if ~isfield(gdat_data.gdat_params,'time_interval')
gdat_data.gdat_params.time_interval = [];
end
- [aa,bb]=unix(['ssh ' '-o "StrictHostKeyChecking no" sxd3d20.d3d.ipp.mpg.de WhichSX ' num2str(shot) ' ' ...
- 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';
+ exp_name_eff = 'D3D';
icount = 0;
nnth = 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')
gdat_data.gdat_params.freq = nnth;
end
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
- if ~isempty(chords_ok_diags(i).diag) && ~isempty(chords_ok_diags(i).chord)
- [a,e]=rdaD3D_eff(shot,chords_ok_diags(i).diag,chords_ok_diags(i).chord,exp_name_eff,gdat_data.gdat_params.time_interval);
- else
- a.data = [];
- a.t = [];
- end
+ tree = 'spectroscopy';
+ ichord = gdat_data.gdat_params.camera(i);
+ diagname = ['sxr:' gdat_data.gdat_params.source ':' gdat_data.gdat_params.source num2str(ichord,'%.2d')];
+ a = gdat_d3d(shot,{tree,diagname});
if ~isempty(a.data)
icount = icount + 1;
if icount == 1
% 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.units = a.units;
- if strcmp(dim1_len,'from_chord_nb')
- gdat_data.data = NaN*ones(max(chords_ok_nb),length(gdat_data.t));
- gdat_data.dim{1} = [1:max(chords_ok_nb)]; % simpler to have index corresponding to chord number, except for central
- else
- gdat_data.data = NaN*ones(length(chords_ok_nb),length(gdat_data.t));
- gdat_data.dim{1} = chords_ok_nb;
- end
+ gdat_data.data = NaN*ones(max(gdat_data.gdat_params.camera),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} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
gdat_data.dimunits = [{'chord nb'}; {'s'}];
- gdat_data.data_fullpath = ['sxr from source=''' gdat_data.gdat_params.source ' uses shotfiles SX' ...
- setdiff(unique(strvcat(chords_ok_diags(:).diag)),'SX')'];
+ gdat_data.data_fullpath = ['sxr from source=''' gdat_data.gdat_params.source ''''];
gdat_data.label = ['SXR/' upper(gdat_data.gdat_params.source)];
end
try
- if strcmp(dim1_len,'from_chord_nb')
- gdat_data.data(chords_ok_nb(i),:) = a.data(1:nnth:end);
- else
- gdat_data.data(i,:) = a.data(1:nnth:end);
- end
+ gdat_data.data(ichord,:) = a.data(1:nnth:end);
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
else
% add fields not yet defined in case all cases have empty data
end
end
- gdat_data.chords = chords_ok_nb;
- if length(gdat_data.dim)>=1
- gdat_data.x = gdat_data.dim{1};
- end
+ gdat_data.chords = gdat_data.gdat_params.camera;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'transp'}
@@ -1936,36 +1644,3 @@ return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);
diff --git a/crpptbx/D3D/get_eqdsk_d3d.m b/crpptbx/D3D/get_eqdsk_d3d.m
new file mode 100644
index 0000000000000000000000000000000000000000..5fee893e411e5775a2381594e9ba6019a61029c3
--- /dev/null
+++ b/crpptbx/D3D/get_eqdsk_d3d.m
@@ -0,0 +1,128 @@
+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;
diff --git a/crpptbx/D3D/get_signal_d3d.m b/crpptbx/D3D/get_signal_d3d.m
index e1e1a7ebc13edb3d9c5351f1fa26a57d56271ddf..e0d0a0d475e9da566d5a3ac5a1b34f2870a1ced1 100644
--- a/crpptbx/D3D/get_signal_d3d.m
+++ b/crpptbx/D3D/get_signal_d3d.m
@@ -3,26 +3,27 @@ function [gdat_data,error_status] = get_signal_d3d(signal,gdat_data);
gdat_data.data = mdsvalue(signal);
error_status = 1;
+gdat_data.dim = [];
+gdat_data.dimunits = [];
+gdat_data.t = [];
+gdat_data.x = [];
nbdims=numel(setdiff(size(gdat_data.data),1));
if ~isempty(gdat_data.data) && ~ischar(gdat_data.data)
gdat_data.units = mdsvalue(['units_of(' signal ')']);
- for i=1:nbdims
- gdat_data.dim{i} = mdsvalue(['dim_of(' signal ',' num2str(i-1) ')']);
- % gdat_data.dimunits{i} = mdsvalue(['dimunits_of(' signal ',' num2str(i-1) ')']);
- gdat_data.dimunits{i} = 'to get';
+ if prod(size(gdat_data.data)) > 1
+ for i=1:nbdims
+ gdat_data.dim{i} = mdsvalue(['dim_of(' signal ',' num2str(i-1) ')']);
+ % 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
- 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