diff --git a/crpptbx/AUG/loadAUGdata.m b/crpptbx/AUG/loadAUGdata.m
index 038bdc2f6924e452f83dcbf18770bbed974db47f..897ac2802e0abb6c2bcd0c911e6c9d329ef0bd71 100644
--- a/crpptbx/AUG/loadAUGdata.m
+++ b/crpptbx/AUG/loadAUGdata.m
@@ -1081,54 +1081,68 @@ switch AUGkeywrdcase{index}
[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))));
- 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(equil.qvalue(:,it)~=0);
- qeff=equil.qvalue(ijk(end:-1:1),it);
- rhoeff=equil.rhopolnorm(3:2:2*length(ijk)+1,it);
- if length(ijk)>2
- ij_nonan=find(~isnan(equil.rhopolnorm(:,it)));
- qfit = zeros(size(equil.rhopolnorm(:,it)));
- qfit(ij_nonan)=interpos([0.; rhoeff],[qeff(1) ;qeff],equil.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff))]);
+ % check NTIME ok with qvalue
+ if NTIME== size(qpsi.value,2)
+ 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
- qfit = zeros(size(equil.rhopolnorm(:,it)));
+ equil.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
end
- equil.qvalue(:,it) = qfit;
- end
- % 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);
- %
+ % 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))));
+ 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(equil.qvalue(:,it)~=0);
+ qeff=equil.qvalue(ijk(end:-1:1),it);
+ rhoeff=equil.rhopolnorm(3:2:2*length(ijk)+1,it);
+ if length(ijk)>2
+ ij_nonan=find(~isnan(equil.rhopolnorm(:,it)));
+ qfit = zeros(size(equil.rhopolnorm(:,it)));
+ qfit(ij_nonan)=interpos([0.; rhoeff],[qeff(1) ;qeff],equil.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff))]);
+ else
+ qfit = zeros(size(equil.rhopolnorm(:,it)));
+ end
+ equil.qvalue(:,it) = qfit;
+ end
+ % 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
+ else
+ disp('***********************************************************************')
+ disp('***********************************************************************')
+ disp(['problem with nb of time points: NTIME for parameter is: ' num2str(NTIME) ' while qpsi dim= ' num2str(size(qpsi.value))])
+ disp('***********************************************************************')
+ disp('***********************************************************************')
+ equil.qvalue = [];
+ equil.qpsi=qpsi;
+ equil.PFL=psi_tree;
+ equil.TFLx=phi_tree;
+ equil.rhopolnorm=[];
end
equil.x = equil.rhopolnorm;
%