From e868989f7d6481b59d939bbce5f4f8fa51ca777f Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Wed, 30 Apr 2014 14:14:57 +0000
Subject: [PATCH] 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 | 85 +++++++++++++++++++++++++++++++++++++++++
1 file changed, 85 insertions(+)
create mode 100644 crpptbx/addpsi_fromRZ.m
diff --git a/crpptbx/addpsi_fromRZ.m b/crpptbx/addpsi_fromRZ.m
new file mode 100644
index 00000000..d3aaa784
--- /dev/null
+++ b/crpptbx/addpsi_fromRZ.m
@@ -0,0 +1,85 @@
+
+% $$$ 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
--
GitLab