add R and rho in ece, eced, ece_rho and eced_rho gdat calls

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@4449 d63d8f72-b253-0410-a779-e742ad2e26cf

crpptbx/AUG/loadAUGdata.m

@@ -83,7 +83,9 @@ function [trace,error,varargout]=loadAUGdata(shot,data_type,varargin)
 % sxb : 'SXB/J' chords
 % sxf : 'SXF/I' chords
 % ece : 
+% ece_rho : 
 % eced : 
+% eced_rho : 
 % eced_rmd : 
 % Halpha : 
 % pgyro : for each gyrotrons, power, freq, etc (ask for more)
@@ -185,9 +187,15 @@ if size(data_type_eff,1)==1
   if ~isempty(strmatch(data_type_eff_noext,[{'ECE'}],'exact'))
     data_type_eff_noext='ece';
   end
+  if ~isempty(strmatch(lower(data_type_eff_noext),[{'ece_rho'}],'exact'))
+    data_type_eff_noext='ece_rho';
+  end
   if ~isempty(strmatch(upper(data_type_eff_noext),[{'ECED'}],'exact'))
     data_type_eff_noext='eced';
   end
+  if ~isempty(strmatch(upper(data_type_eff_noext),[{'ECED_RHO'}],'exact'))
+    data_type_eff_noext='eced_rho';
+  end
   if ~isempty(strmatch(upper(data_type_eff_noext),[{'ECED_RMD'}],'exact'))
     data_type_eff_noext='eced_rmd';
   end
@@ -261,7 +269,8 @@ end
 AUGkeywrdall=[{'Ip'} {'b0'} {'zmag'} {'rmag'}  {'rgeo'} {'zgeo'} {'rcont'} {'zcont'} {'vol'} {'qrho'} {'qrho_fpp'} {'q0'} {'q95'} {'kappa'} ...
 	      {'delta'} {'deltatop'} {'deltabot'} {'neint'} {'ne'} {'te'}  ...
 	      {'nerho'} {'neterho'} {'terho'} {'cxrs'} {'cxrs_rho'} {'equil'} {'equil_fpp'} {'equil_eqm'} ...
-	      {'equil_eqr'} {'sxr'} {'sxR'} {'sxb'} {'sxf'} {'ssx_g'} {'ssx_h'} {'ssx_i'} {'ssx_j'} {'ssx'} {'ece'} {'eced'} {'eced_rmd'} {'Halpha'} {'pgyro'} {'powers'}];
+	      {'equil_eqr'} {'sxr'} {'sxR'} {'sxb'} {'sxf'} {'ssx_g'} {'ssx_h'} {'ssx_i'} {'ssx_j'} {'ssx'} ...
+	      {'ece'} {'ece_rho'} {'eced'} {'eced_rho'} {'eced_rmd'} {'Halpha'} {'pgyro'} {'powers'}];
 AUGsig.iip=strmatch('Ip',AUGkeywrdall,'exact');
 AUGsig.ib0=strmatch('b0',AUGkeywrdall,'exact');
 AUGsig.izmag=strmatch('zmag',AUGkeywrdall,'exact');
@@ -302,6 +311,8 @@ AUGsig.issx_j=strmatch('ssx_j',AUGkeywrdall,'exact');
 AUGsig.issx=strmatch('ssx',AUGkeywrdall,'exact');
 AUGsig.iece=strmatch('ece',AUGkeywrdall,'exact');
 AUGsig.ieced=strmatch('eced',AUGkeywrdall,'exact');
+AUGsig.iece_rho=strmatch('ece_rho',AUGkeywrdall,'exact');
+AUGsig.ieced_rho=strmatch('eced_rho',AUGkeywrdall,'exact');
 AUGsig.ieced_rmd=strmatch('eced_rmd',AUGkeywrdall,'exact');
 AUGsig.iHalpha=strmatch('Halpha',AUGkeywrdall,'exact');
 AUGsig.ipgyro=strmatch('pgyro',AUGkeywrdall,'exact');
@@ -334,6 +345,8 @@ AUGkeywrdcase(AUGsig.issx_j)=AUGkeywrdall(AUGsig.issx_j);
 AUGkeywrdcase(AUGsig.issx)=AUGkeywrdall(AUGsig.issx);
 AUGkeywrdcase(AUGsig.iece)=AUGkeywrdall(AUGsig.iece);
 AUGkeywrdcase(AUGsig.ieced)=AUGkeywrdall(AUGsig.ieced);
+AUGkeywrdcase(AUGsig.iece_rho)=AUGkeywrdall(AUGsig.iece_rho);
+AUGkeywrdcase(AUGsig.ieced_rho)=AUGkeywrdall(AUGsig.ieced_rho);
 AUGkeywrdcase(AUGsig.icxrs)=AUGkeywrdall(AUGsig.icxrs);
 AUGkeywrdcase(AUGsig.icxrs_rho)=AUGkeywrdall(AUGsig.icxrs_rho);
 AUGkeywrdcase(AUGsig.ipgyro)=AUGkeywrdall(AUGsig.ipgyro); % idem
@@ -364,6 +377,8 @@ AUGsiglocation(:,AUGsig.ideltabot)={''; ''};
 AUGsiglocation(:,AUGsig.ineint)={'DCN'; 'H-1'};
 AUGsiglocation(:,AUGsig.iece)={'CEC'; 'Trad-A'};
 AUGsiglocation(:,AUGsig.ieced)={'CEC'; 'Trad-A'}; % ECED
+AUGsiglocation(:,AUGsig.iece_rho)={'CEC'; 'Trad-A'};
+AUGsiglocation(:,AUGsig.ieced_rho)={'CEC'; 'Trad-A'}; % ECED
 AUGsiglocation(:,AUGsig.ieced_rmd)={'RMD'; 'Trad-A'}; % ECED_RMD
 AUGsiglocation(:,AUGsig.iHalpha)={'POT'; 'ELMa-Han'};
 
@@ -373,6 +388,7 @@ AUGexplocation=cell(size(AUGkeywrdall,2),1);
 AUGexplocation(:)={'AUGD'};
 % cases with a "private" shotfile:
 AUGexplocation(AUGsig.ieced)={'ECED'};
+AUGexplocation(AUGsig.ieced_rho)={'ECED'};
 AUGexplocation(AUGsig.ieced_rmd)={'ECED'};
 
 % initialize order of substructures and allows just a "return" if data empty
@@ -658,7 +674,7 @@ switch AUGkeywrdcase{index}
     end
     
   %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-  case {'ece','eced'}
+  case {'ece','eced','ece_rho','eced_rho'}
     %  LOAD MULTI CHANNEL DATA    
     %  load AUG ece data
 
@@ -705,21 +721,65 @@ switch AUGkeywrdcase{index}
     trace.dim{2}=trace.t;
     % get R
     [aR,e]=rdaAUG_eff(shot,ppftype,'R-A',shotfile_exp_eff,timerange);
-    domatchRtime=1;
+    [aZ,e]=rdaAUG_eff(shot,ppftype,'z-A',shotfile_exp_eff,timerange);
+    domatchRtime=0;
     if domatchRtime
       % interpolate R structure on ece data time array, to ease plot vs R
       for i=starti:endi
 	radius.data(i,:) = interp1(aR.t,aR.data(i,:),trace.t);
+	zheight.data(i,:) = interp1(aZ.t,aZ.data(i,:),trace.t);
       end
       radius.t=trace.t;
+      zheight.t=trace.t;
     else
       radius.data = aR.data;
       radius.t=aR.t;
+      zheight.data = aZ.data;
+      zheight.t=aR.t;
     end
     ij=find(trace.data==0);
     trace.data(ij)=NaN;
     varargout{1}={radius};
     trace.R=radius;
+    trace.Z=zheight;
+
+    if strcmp(AUGkeywrdcase{index},'ece_rho') || strcmp(AUGkeywrdcase{index},'eced_rho')
+      equil=gdat(shot,'equil',0);
+      inb_chord_ece=size(trace.R.data,1);
+      inb_time_ece=size(trace.R.data,2);
+      psi_out = NaN*ones(inb_chord_ece,inb_time_ece);
+      rhopsinorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
+      rhotornorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
+      rhovolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
+      % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
+      time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-1)];
+      iok=find(aR.data(:,1)>0);
+      for itequil=1:length(time_equil)-1
+	rr=equil.Rmesh(:,itequil);
+	zz=equil.Zmesh(:,itequil);
+	psirz_in = equil.psi2D(:,:,itequil);
+	it_ece_inequil = find(trace.R.t>=time_equil(itequil) & trace.R.t<=time_equil(itequil+1));
+	if ~isempty(it_ece_inequil)
+	  rout=trace.R.data(iok,it_ece_inequil);
+	  ijok=find(~isnan(rout));
+	  if ~isempty(ijok)
+	    zout=trace.Z.data(iok,it_ece_inequil);
+	    psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
+	    psi_out(iok,it_ece_inequil) = reshape(psi_at_routzout,length(iok),length(it_ece_inequil));
+	    rhopsinorm_out(iok,it_ece_inequil) = sqrt((psi_out(iok,it_ece_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
+	    for it_cx=1:length(it_ece_inequil)
+	      rhotornorm_out(iok,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out(iok,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
+	      rhovolnorm_out(iok,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out(iok,it_ece_inequil(it_cx)),-3,[2 2],[0 1]);
+	    end
+	  end
+	end
+      end
+      trace.rhos.psi_on_rztime = psi_out;
+      trace.rhos.rhopsinorm_on_rztime = rhopsinorm_out;
+      trace.rhos.rhotornorm_on_rztime = rhotornorm_out;
+      trace.rhos.rhovolnorm_on_rztime = rhovolnorm_out;
+      trace.rhos.t = trace.R.t;
+    end
 
   %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
   case {'cxrs', 'cxrs_rho'}