add example to calculate psidiag(Rdiag,Zdiag) at specific diagnostic...

add example to calculate psidiag(Rdiag,Zdiag) at specific diagnostic Rdiag,Zdiag values from equilibrium psi(R,Z), as well as rhopol, rhotor etc. This is used in gdat,AUG now

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

crpptbx/addpsi_fromRZ.m

+
+% $$$ shot=30382;
+% $$$ time=3.;
+% $$$ 
+% $$$ cxrs=gdat(shot,'cxrs',0);
+% $$$ equil=gdat(shot,'equil',0);
+
+if ~exist('equil')
+  disp('no equil')
+  return
+end
+if ~exist('cxrs')
+  disp('no cxrs')
+  return
+end
+
+inb_chord_cxrs=size(cxrs.r_time,1);
+inb_time_cxrs=size(cxrs.r_time,2);
+psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
+rhopsinorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
+rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
+rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
+% 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(cxrs.r_time(:,1)>0);
+tic
+for itequil=1:length(time_equil)-1
+  rr=equil.Rmesh(:,itequil);
+  zz=equil.Zmesh(:,itequil);
+  psirz_in = equil.psi2D(:,:,itequil);
+  it_cxrs_inequil = find(cxrs.t>=time_equil(itequil) & cxrs.t<=time_equil(itequil+1));
+  if ~isempty(it_cxrs_inequil)
+    rout=cxrs.r_time(iok,it_cxrs_inequil);
+    zout=cxrs.z_time(iok,it_cxrs_inequil);
+    psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
+    psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil));
+    rhopsinorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
+    for it_cx=1:length(it_cxrs_inequil)
+      rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
+      rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
+    end
+  end
+end
+toc
+
+cxrs_all = cxrs;
+cxrs_all.psi_on_rztime = psi_out;
+cxrs_all.rhopsinorm_on_rztime = rhopsinorm_out;
+cxrs_all.rhotornorm_on_rztime = rhotornorm_out;
+cxrs_all.rhovolnorm_on_rztime = rhovolnorm_out;
+
+% $$$ iok=find(cxrs.r_time(:,1)>0);
+% $$$ rout=cxrs.r_time(iok,1);
+% $$$ zout=cxrs.z_time(iok,1);
+% $$$ 
+% $$$ % [farray_out,dfarrdx_out,dfarrdy_out,d2farrdx2_out,d2farrdy2_out,d2farrdxdy_out]= ...
+% $$$ [farray_out]= interpos2Dcartesian(rr',zz',farray_in,rout,zout);
+% $$$ 
+% $$$ figure;
+% $$$ contour(rr,zz,farray_in(:,:,1)',100);
+% $$$ hold on
+% $$$ plot(rout,zout,'*')
+% $$$ axis equal
+% $$$ colorbar
+% $$$ 
+% $$$ figure;
+% $$$ contour(rr,zz,farray_in(:,:,30)',linspace(equil.psi_axis(30),0..*equil.psi_axis(30)+1.0.*equil.psi_lcfs(30),100));
+% $$$ hold on
+% $$$ plot(rout,zout,'*')
+% $$$ axis equal
+% $$$ colorbar
+% $$$ 
+% $$$ [zzz itime]=min(abs(equil.t-time-1.5));
+% $$$ [zzz iz]=min(abs(zz-0.05));
+% $$$ figure
+% $$$ plot(rr,farray_in(:,iz,itime),'*')
+% $$$ hold on
+% $$$ plot(rout(iok),farray_out(iok,1,itime),'ro')
+% $$$ 
+% $$$ figure;
+% $$$ contour(rr,zz,farray_in(:,:,itime)',farray_out(iok,1,itime));
+% $$$ hold on
+% $$$ plot(rout,zout,'*')
+% $$$ axis equal
+% $$$ colorbar