From a96c27fe6bb6157cd7c0a4309e3b8e17bd6e54f3 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Tue, 24 Jun 2014 14:22:09 +0000
Subject: [PATCH] fix present problem of EQI, EQR, etc from within loadAUG
(time base not consistent) with a method which should be OK when OK, taking
NTIME first points
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@4468 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/AUG/loadAUGdata.m | 174 +++++++++++++++++++++++---------------
1 file changed, 105 insertions(+), 69 deletions(-)
diff --git a/crpptbx/AUG/loadAUGdata.m b/crpptbx/AUG/loadAUGdata.m
index 680ec9d7..aa85adf3 100644
--- a/crpptbx/AUG/loadAUGdata.m
+++ b/crpptbx/AUG/loadAUGdata.m
@@ -460,10 +460,23 @@ if size(data_type_eff,1)==2
AUGsiglocation(1:2,end+1)=[data_type_eff(1) ; {data_type_eff_noext}];
AUGexplocation{end+1}=shotfile_exp;
AUGsigtimeindx(end+1)=0;
- elseif ~strcmp(AUGkeywrdcase{index},'simplereaddata')
- msgbox(['Problem in loadAUGdata with data_type_eff = ' char(data_type_eff(end)) ...
- '. Full paths of nodes should only be for case simplereaddata'],'in loadAUGdata','error')
- error('in loadAUGdata')
+ AUGkeywrdcase{index}
+ elseif length(index)>1 || ~strcmp(AUGkeywrdcase{index},'simplereaddata')
+% $$$ msgbox(['Problem in loadAUGdata with data_type_eff = ' char(data_type_eff(end)) ...
+% $$$ '. Full paths of nodes should only be for case simplereaddata'],'in loadAUGdata','error')
+% $$$ error('in loadAUGdata')
+ % chose a full path which coincides with a keyword
+ disp('the node path seems to coincide with the following keywords, which might be better to use:')
+ disp(' ')
+ disp(AUGkeywrdcase(index))
+ disp('proceed as is in any case')
+ index=length(AUGkeywrdall)+1;
+ AUGkeywrdall(end+1)={'new'};
+ AUGkeywrdcase(end+1)={'simplereaddata'};
+ AUGsiglocation(1:2,end+1)=[data_type_eff(1) ; {data_type_eff_noext}];
+ AUGexplocation{end+1}=shotfile_exp;
+ AUGsigtimeindx(end+1)=0;
+ AUGkeywrdcase{index}
end
else
index=strmatch(data_type_eff_noext,AUGkeywrdall,'exact');
@@ -1059,6 +1072,9 @@ switch AUGkeywrdcase{index}
case {'equil', 'equil_fpp', 'equil_eqm', 'equil_eqr', 'equil_eqh', 'qrho', 'qrho_fpp'}
+ % In shotfiles time is 1st index in 2D arrays while gdat convention is to have time as last dimension
+ % so change equil.xxx to (rho,time) on the way...
+
shotfile_exp_eff = AUGexplocation{index};
if strcmp(AUGkeywrdcase{index},'equil_fpp') || strcmp(AUGkeywrdcase{index},'qrho_fpp')
@@ -1079,90 +1095,110 @@ switch AUGkeywrdcase{index}
NTIME_par = sf2par(DIAG,shot,'NTIME','PARMV');
NTIME = NTIME_par.value; % nb of points
Lpf_par = rdaAUG_eff(shot,DIAG,'Lpf',shotfile_exp_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);
- Lpf_tot = Lpf_par.value; % nb of points: 100000*nb_SOL points + nb_core
+ if (NTIME < NTIME_Lpf)
+ disp('WARNING: nb of times points smaller then equil results, use first NTIME points')
+ elseif (NTIME > NTIME_Lpf)
+ disp('ERROR: nb of times points LARGER then equil results')
+ disp('this is unexpected, so stop there and ask Olivier.Sauter@epfl.ch')
+ 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
-% since Lpf depends on time, need to load all first and then loop over time for easier mapping
+ % 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);
+ 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
+ itotransposeback = 0;
+ end
+ qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback);
[psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',shotfile_exp);
+ psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
[phi_tree,e]=rdaAUG_eff(shot,DIAG,'TFLx',shotfile_exp);
+ phi_tree=adapt_rda(phi_tree,NTIME,ndimrho,itotransposeback);
[Vol,e]=rdaAUG_eff(shot,DIAG,'Vol',shotfile_exp);
+ Vol=adapt_rda(Vol,NTIME,2*ndimrho,itotransposeback);
[Area,e]=rdaAUG_eff(shot,DIAG,'Area',shotfile_exp);
+ Area=adapt_rda(Area,NTIME,2*ndimrho,itotransposeback);
[Ri,e]=rdaAUG_eff(shot,DIAG,'Ri',shotfile_exp);
+ Ri=adapt_rda(Ri,NTIME,M_Rmesh,itotransposeback);
[Zj,e]=rdaAUG_eff(shot,DIAG,'Zj',shotfile_exp);
+ Zj=adapt_rda(Zj,NTIME,N_Zmesh,itotransposeback);
[PFM_tree,e]=rdaAUG_eff(shot,DIAG,'PFM',shotfile_exp);
+ PFM_tree=adaptPFM_rda(PFM_tree,M_Rmesh,N_Zmesh,NTIME);
[Pres,e]=rdaAUG_eff(shot,DIAG,'Pres',shotfile_exp);
+ Pres=adapt_rda(Pres,NTIME,2*ndimrho,itotransposeback);
[FFP,e]=rdaAUG_eff(shot,DIAG,'FFP',shotfile_exp);
- if isempty(FFP.value); FFP.value=NaN*ones(NTIME,max(Lpf1_t)); end
+ if ~isempty(FFP.value)
+ FFP=adapt_rda(FFP,NTIME,ndimrho,itotransposeback);
+ else
+ FFP.value=NaN*ones(NTIME,max(Lpf1_t));
+ end
[LPFx,e]=rdaAUG_eff(shot,DIAG,'LPFx',shotfile_exp);
+ LPFx.value=LPFx.value(1:NTIME); LPFx.data=LPFx.value; LPFx.t=LPFx.t(1:NTIME);
[PFxx,e]=rdaAUG_eff(shot,DIAG,'PFxx',shotfile_exp);
+ PFxx=adapt_rda(PFxx,NTIME,max(LPFx.value)+1,itotransposeback);
[RPFx,e]=rdaAUG_eff(shot,DIAG,'RPFx',shotfile_exp);
+ RPFx=adapt_rda(RPFx,NTIME,max(LPFx.value)+1,itotransposeback);
[zPFx,e]=rdaAUG_eff(shot,DIAG,'zPFx',shotfile_exp);
+ 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 = 40; % Note could just put a NaN on LCFS value since ill-defined when diverted
- % 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);
+ max_qValue = 40.0; % 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))]);
else
- equil.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
+ qfit = zeros(size(equil.rhopolnorm(:,it)));
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))]);
- 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=[];
+ 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
equil.x = equil.rhopolnorm;
%
@@ -1175,7 +1211,7 @@ switch AUGkeywrdcase{index}
equil.t = equil_time.value(1:NTIME);
%
[IpiPSI,e]=rdaAUG_eff(shot,DIAG,'IpiPSI',shotfile_exp);
- equil.Ip = IpiPSI.value';
+ equil.Ip = IpiPSI.value(1:NTIME);
%
equil.data = equil.qvalue; % put q in data
equil.unit=[]; % not applicable
--
GitLab