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

add protection when reading EQxx with bad time arrays, returns warning and qpsi PFL TFLx as is 

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@4460 d63d8f72-b253-0410-a779-e742ad2e26cf
parent 50bfad9e
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
crpptbx/AUG/loadAUGdata.m
+ 60
− 46
View file @ 2929c443
......@@ -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;
%
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment