diff --git a/crpptbx/AUG/loadAUGdata.m b/crpptbx/AUG/loadAUGdata.m
index aa85adf3a74c648e36208f55a585cdb45a84150a..8a4f122dcd158d0e6fdf94cccaeb3c338d927e23 100644
--- a/crpptbx/AUG/loadAUGdata.m
+++ b/crpptbx/AUG/loadAUGdata.m
@@ -16,6 +16,7 @@ function [trace,error,varargout]=loadAUGdata(shot,data_type,varargin)
 %     gdat(15133,'AUGD/MAG/Ipa',1,'AUG') % to specify experiment explicitely like in:
 %     gdat(30230,'ECED/RMD/Trad-A',1,'AUG') % 
 %     gdat(30230,'ECED/CEC/Trad-A',1,'AUG') % 
+%     gdat(31053,'MSP/IGP_07Co',1); % N gas opening at -6s (03 as well)
 %
 % INPUT:
 % shot: shot number
@@ -729,7 +730,6 @@ switch AUGkeywrdcase{index}
     %  load AUG ece data
 
     shotfile_exp_eff = AUGexplocation{index};
-
     if  nargin>=3 & ~isempty(varargin{1})
       starti=varargin{1}(1);
       endi=varargin{1}(2);
@@ -770,8 +770,16 @@ switch AUGkeywrdcase{index}
     trace.data=trace.data(:,1:nth:end);
     trace.dim{2}=trace.t;
     % get R
-    [aR,e]=rdaAUG_eff(shot,ppftype,'R-A',shotfile_exp_eff,timerange);
-    [aZ,e]=rdaAUG_eff(shot,ppftype,'z-A',shotfile_exp_eff,timerange);
+    try
+      [aR,e]=rdaAUG_eff(shot,ppftype,'R-A',shotfile_exp_eff,timerange);
+    catch
+    end
+    try
+      % [aZ,e]=rdaAUG_eff(shot,ppftype,'z-A',shotfile_exp_eff,timerange); (problem now)
+      aZ=aR;
+      aZ.data=ones(size(aR.data));
+    catch
+    end
     domatchRtime=0;
     if domatchRtime
       % interpolate R structure on ece data time array, to ease plot vs R
@@ -1119,6 +1127,8 @@ switch AUGkeywrdcase{index}
       itotransposeback = 0;
     end
     qpsi=adapt_rda(qpsi,NTIME,ndimrho,itotransposeback);
+    ijnan=find(isnan(qpsi.value));
+    qpsi.value(ijnan)=0;
     [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);
@@ -1135,12 +1145,28 @@ switch AUGkeywrdcase{index}
     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);
+    [Jpol,e]=rdaAUG_eff(shot,DIAG,'Jpol',shotfile_exp);
+    Jpol=adapt_rda(Jpol,NTIME,2*ndimrho,itotransposeback);
     [FFP,e]=rdaAUG_eff(shot,DIAG,'FFP',shotfile_exp);
     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')
+      [Rinv,e]=rdaAUG_eff(shot,DIAG,'Rinv',shotfile_exp);
+      Rinv=adapt_rda(Rinv,NTIME,ndimrho,itotransposeback);
+      [R2inv,e]=rdaAUG_eff(shot,DIAG,'R2inv',shotfile_exp);
+      R2inv=adapt_rda(R2inv,NTIME,ndimrho,itotransposeback);
+      [Bave,e]=rdaAUG_eff(shot,DIAG,'Bave',shotfile_exp);
+      Bave=adapt_rda(Bave,NTIME,ndimrho,itotransposeback);
+      [B2ave,e]=rdaAUG_eff(shot,DIAG,'B2ave',shotfile_exp);
+      B2ave=adapt_rda(B2ave,NTIME,ndimrho,itotransposeback);
+      [FTRA,e]=rdaAUG_eff(shot,DIAG,'FTRA',shotfile_exp);
+      FTRA=adapt_rda(FTRA,NTIME,ndimrho,itotransposeback);
+    else
+      Rinv.value=[]; R2inv.value=[]; Bave.value=[]; B2ave.value=[]; FTRA.value=[];
+    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);
@@ -1155,7 +1181,9 @@ switch AUGkeywrdcase{index}
       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
+      ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially
+      % set NaNs to zeroes
+      if qpsi.value(ijok(1),1)<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);
@@ -1169,12 +1197,17 @@ switch AUGkeywrdcase{index}
 	% 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
+	  % now shots have non-zero axis values in eqr
+	  rhoeff=equil.rhopolnorm(ijk,it);
+	  qeff=equil.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(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))]);
+	  qfit(ij_nonan)=interpos(rhoeff,qeff,equil.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]);
 	else
 	  qfit = zeros(size(equil.rhopolnorm(:,it)));
 	end
@@ -1194,10 +1227,42 @@ switch AUGkeywrdcase{index}
       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
+      if ~isempty(Jpol.value)
+	equil.jpol(:,it)=Jpol.value(it,2*Lpf1-1:-2:1)'; 
+	equil.djpolpsi(:,it)=Jpol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi
+      else
+	equil.jpol = [];
+	equil.djpolpsi = [];
+      end
       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);
+      if ~isempty(Rinv.value)
+	equil.rinv(:,it) = Rinv.value(it,Lpf1:-1:1)';
+      else
+	equil.rinv = [];
+      end
+      if ~isempty(R2inv.value)
+	equil.r2inv(:,it) = R2inv.value(it,Lpf1:-1:1)';
+      else
+	equil.r2inv = [];
+      end
+      if ~isempty(Bave.value)
+	equil.bave(:,it) = Bave.value(it,Lpf1:-1:1)';
+      else
+	equil.bave = [];
+      end
+      if ~isempty(B2ave.value)
+	equil.b2ave(:,it) = B2ave.value(it,Lpf1:-1:1)';
+      else
+	equil.b2ave = [];
+      end
+      if ~isempty(FTRA.value)
+	equil.ftra(:,it) = FTRA.value(it,Lpf1:-1:1)';
+      else
+	equil.ftra = [];
+      end
       %
     end
     equil.x = equil.rhopolnorm;