diff --git a/crpptbx/AUG/loadAUGdata.m b/crpptbx/AUG/loadAUGdata.m
index 5eafeb6680d4729adf265eeb17ad0240279df1a0..9e535de1c6848684e808d1d5b045d680512ae3e9 100644
--- a/crpptbx/AUG/loadAUGdata.m
+++ b/crpptbx/AUG/loadAUGdata.m
@@ -144,9 +144,9 @@ end
% all keywords and corresponding case to run below
AUGkeywrdall=[{'Ip'} {'zmag'} {'rmag'} {'rcont'} {'zcont'} {'vol'} {'qrho'} {'qrho_cliste'} {'q95'} {'kappa'} ...
- {'delta'} {'deltatop'} {'deltabot'} {'neint'} ...
- {'ne'} {'te'} {'nerho'} {'terho'} ...
- {'sxr'} {'sxR'} {'sxb'} {'ece'} {'eced'} {'Halpha'}];
+ {'delta'} {'deltatop'} {'deltabot'} {'neint'} ...
+ {'ne'} {'te'} {'nerho'} {'terho'} {'equil'} {'equil_cliste'} ...
+ {'sxr'} {'sxR'} {'sxb'} {'ece'} {'eced'} {'Halpha'}];
AUGsig.iip=strmatch('Ip',AUGkeywrdall,'exact');
AUGsig.izmag=strmatch('zmag',AUGkeywrdall,'exact');
AUGsig.irmag=strmatch('rmag',AUGkeywrdall,'exact');
@@ -155,6 +155,8 @@ AUGsig.izcont=strmatch('zcont',AUGkeywrdall,'exact');
AUGsig.ivol=strmatch('vol',AUGkeywrdall,'exact');
AUGsig.iqrho=strmatch('qrho',AUGkeywrdall,'exact');
AUGsig.iqrho_cliste=strmatch('qrho_cliste',AUGkeywrdall,'exact');
+AUGsig.iequil=strmatch('equil',AUGkeywrdall,'exact');
+AUGsig.iequil_cliste=strmatch('equil_cliste',AUGkeywrdall,'exact');
AUGsig.iq95=strmatch('q95',AUGkeywrdall,'exact');
AUGsig.ikappa=strmatch('kappa',AUGkeywrdall,'exact');
AUGsig.idelta=strmatch('delta',AUGkeywrdall,'exact');
@@ -178,6 +180,8 @@ AUGkeywrdcase=cell(size(AUGkeywrdall));
AUGkeywrdcase(:)={'simplereaddata'};
AUGkeywrdcase(AUGsig.iqrho)=AUGkeywrdall(AUGsig.iqrho); % special as efit q on psi
AUGkeywrdcase(AUGsig.iqrho_cliste)=AUGkeywrdall(AUGsig.iqrho_cliste); % special as efit q on psi
+AUGkeywrdcase(AUGsig.iequil)=AUGkeywrdall(AUGsig.iequil); % special as efit q on psi
+AUGkeywrdcase(AUGsig.iequil_cliste)=AUGkeywrdall(AUGsig.iequil_cliste); % special as efit q on psi
%AUGkeywrdcase(AUGsig.idelta)=AUGkeywrdall(AUGsig.idelta); % special as average of triu and tril
AUGkeywrdcase(AUGsig.ine)=AUGkeywrdall(AUGsig.ine); % special as adds error bars
AUGkeywrdcase(AUGsig.ite)=AUGkeywrdall(AUGsig.ite); % idem
@@ -486,47 +490,91 @@ switch AUGkeywrdcase{index}
trace.dimunits{1} = a.area.name ;
trace.dimunits{2} = a.time_aug.unit;
- case {'qrho', 'qrho_cliste'}
+ case {'equil', 'equil_cliste', 'qrho', 'qrho_cliste'}
shotfile_exp_eff = AUGexplocation{index};
- if strcmp(AUGkeywrdcase{index},'qrho')
+ if strcmp(AUGkeywrdcase{index},'equil') || strcmp(AUGkeywrdcase{index},'qrho')
DIAG = 'FPP';
else
DIAG = 'EQI';
end
- [a,e]=rdaAUG_eff(shot,DIAG,'Qpsi',shotfile_exp_eff);
- % Qpsi has inverted channel/time from CEC
- a.value = a.value(:,end:-1:1)';
- a.data = a.value;
- % get x values
- [psi,e]=rdaAUG_eff(shot,DIAG,'PFL',shotfile_exp);
- psi.value=psi.value(:,end:-1:1)';
- psi_axis= (ones(size(psi.value,1),1) * psi.value(1,:));
- a.x = sqrt(abs((psi.value-psi_axis) ./(zeros(size(psi.value))-psi_axis) )); % psi_edge=0 assumed
- a.psi = psi.value;
- % get time values
- [psi_time,e]=rdaAUG_eff(shot,DIAG,'time',shotfile_exp);
- if length(psi_time.value) > length(a.x(1,:))
- a.t = psi_time.value(1:length(a.x(1,:))); % problem with times??
- elseif length(psi_time.value) < length(a.x(1,:))
- a.t = psi_time.value; % problem with times??
- a.t(end+1:length(a.x(1,:))) = psi_time.value(end) + [1:length(a.x(1,:))-length(psi_time.value)].*diff(psi_time.value(end-1:end));
- else
- a.t = psi_time.value;
- end
- % get rhotor values
- [phi,e]=rdaAUG_eff(shot,DIAG,'TFLx',shotfile_exp);
- phi.value=phi.value(:,end:-1:1)';
- a.rhotor = sqrt(abs(phi.value ./ (ones(size(phi.value,1),1) * phi.value(end,:))));
- a.torflux = phi.value;
- % get rhovol values
+ 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',shotfile_exp_eff);
+ Lpf_tot = Lpf_par.value; % nb of points: 100000*nb_SOL points + nb_core
+ Lpf_SOL = fix(Lpf_tot(1:NTIME)/100000);
+ Lpf1_t = mod(Lpf_tot(1:NTIME),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',shotfile_exp_eff);
+ [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',shotfile_exp);
+ [phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',shotfile_exp);
[Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',shotfile_exp);
- Vol.value=Vol.value(:,end-1:-2:1)'; % 2nd index are dV/dpsi
- a.rhovol = sqrt(abs(Vol.value ./ (ones(size(Vol.value,1),1) * Vol.value(end,:))));
- a.volume = Vol.value;
+ [Area,e]=rdaAUG_eff(shot,DIAG,'Area',shotfile_exp);
+ [Ri,e]=rdaAUG_eff(shot,DIAG,'Ri',shotfile_exp);
+ [Zj,e]=rdaAUG_eff(shot,DIAG,'Zj',shotfile_exp);
+ [PFM_tree,e]=rdaAUG_eff(shot,DIAG,'PFM',shotfile_exp);
+ [Pres,e]=rdaAUG_eff(shot,DIAG,'Pres',shotfile_exp);
+ [FFP,e]=rdaAUG_eff(shot,DIAG,'FFP',shotfile_exp);
+ if isempty(FFP.value); FFP.value=NaN*ones(NTIME,max(Lpf1_t)); end
+ [LPFx,e]=rdaAUG_eff(shot,DIAG,'LPFx',shotfile_exp);
+ [PFxx,e]=rdaAUG_eff(shot,DIAG,'PFxx',shotfile_exp);
+ [RPFx,e]=rdaAUG_eff(shot,DIAG,'RPFx',shotfile_exp);
+ [zPFx,e]=rdaAUG_eff(shot,DIAG,'zPFx',shotfile_exp);
+
+ % seems "LCFS" q-value is far too large, limit to some max (when diverted)
+ max_qValue = 40; % 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
+ if qpsi.value(1,2)<0
+ equil.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
+ else
+ equil.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
+ end
+ % get x values
+ equil.psi(:,it)=psi_tree.value(it,Lpf1:-1:1)';
+ equil.psi_axis(it)= equil.psi(1,it);
+ equil.psi_lcfs(it)= equil.psi(end,it);
+ equil.rhopolnorm(:,it) = sqrt(abs((equil.psi(:,it)-equil.psi_axis(it)) ./(equil.psi_lcfs(it)-equil.psi_axis(it))));
+ % get rhotor values
+ equil.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)';
+ equil.rhotornorm(:,it) = sqrt(abs(equil.phi(:,it) ./ equil.phi(end,it)));
+ % get rhovol values
+ equil.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)';
+ equil.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
+ equil.rhovolnorm(:,it) = sqrt(abs(equil.vol(:,it) ./ equil.vol(end,it)));
+ equil.area(:,it)=Area.value(it,2*Lpf1-1:-2:1)';
+ equil.dareadpsi(:,it)=Area.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
+ equil.Rmesh(:,it) = Ri.value(it,1:M_Rmesh);
+ equil.Zmesh(:,it) = Zj.value(it,1:N_Zmesh);
+ equil.psi2D(1:M_Rmesh,1:N_Zmesh,it) = PFM_tree.value(1:M_Rmesh,1:N_Zmesh,it);
+ equil.pressure(:,it)=Pres.value(it,2*Lpf1-1:-2:1)';
+ equil.dpressuredpsi(:,it)=Pres.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
+ equil.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)';
+ equil.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1);
+ equil.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1);
+ equil.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1);
+ %
+ end
+ equil.x = equil.rhopolnorm;
%
- trace = a;
+ [equil_time,e]=rdaAUG_eff(shot,DIAG,'time',shotfile_exp);
+ % get time values
+ equil.t = equil_time.value(1:NTIME);
+ %
+ [IpiPSI,e]=rdaAUG_eff(shot,DIAG,'IpiPSI',shotfile_exp);
+ equil.Ip = IpiPSI.value';
+ %
+ equil.data = equil.qvalue; % put q in data
+ equil.unit=[]; % not applicable
+ equil.name_gdat = [AUGkeywrdcase{index} ' using ' DIAG];
+ trace = equil;
trace.dim{1} = trace.x;
trace.dim{2} = trace.t;
trace.units = trace.unit;
diff --git a/crpptbx/AUG/rdaAUG_eff.m b/crpptbx/AUG/rdaAUG_eff.m
index f8ef57eac2d0191a4da9bba4b112443608f1ffde..5ac675f3ab13d35889d28c20b713aa268b4c1ac4 100644
--- a/crpptbx/AUG/rdaAUG_eff.m
+++ b/crpptbx/AUG/rdaAUG_eff.m
@@ -138,7 +138,7 @@ else
adata.t=[1:size(adata.value,2)];
end
else
- adata.value = adata.value';
+ if length(size(adata.value))<=2; adata.value = adata.value'; end % cannot transpose Nd>2 matrix
if ~isempty(adata.time_aug)
adata.x=[1:prod(size(adata.value))/length(adata_time.value)];
adata.t=adata.time_aug.value;