From 1425bd27d3cf353ef2e424eb991b68f20bcde5e7 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Fri, 1 Feb 2002 13:47:13 +0000
Subject: [PATCH] to calculate a(psi), etc
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@1861 d63d8f72-b253-0410-a779-e742ad2e26cf
---
JET/Ra_rho_t.m | 80 ++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 80 insertions(+)
create mode 100644 JET/Ra_rho_t.m
diff --git a/JET/Ra_rho_t.m b/JET/Ra_rho_t.m
new file mode 100644
index 00000000..f7110235
--- /dev/null
+++ b/JET/Ra_rho_t.m
@@ -0,0 +1,80 @@
+function [R0zmag,azmag,Rminzmag,Rmaxzmag]=Ra_rho_t(sspr,sspi,npts);
+%
+% Rmax(psi^.5),Rmin(psi^.5) and a(psi^.5) on z=zmag on npts equidistant sqrt(psin)
+%
+% Assumes all timing the same: sspr.t, sspi.t, zmag.t, rmag.t
+%
+% cannot use zmag as not good enough with respect to psirz. Better to use min of psirz contours
+%
+% precision within 1e-3 - 1e-4 (as usually min(psi)~a few 1e-4)
+%
+% Note: this also gives Rin and Rout (psi) directly from Rmin and Rmax, or from R0+a or R0-a of rho^2
+%
+
+rhopsi=linspace(0,1,npts)';
+figure(92);clf; % used as dummy for contours calculation
+set(92,'pos',[600 500 200 200])
+set(92,'MenuBar','none');
+hold on
+azmag=sspr; % to get time array and structure same
+azmag.data=zeros(length(rhopsi),length(sspr.t));
+azmag.x=rhopsi;
+azmag.dim{1}=rhopsi;
+azmag.dimunits{1}='rho=sqrt(psin)';
+azmag.name='a(rho,t)=mean minor radius from psi(R,Z) and Ra_rho_t function';
+R0zmag=azmag;
+R0zmag.name='R0(rho,t)=mean major radius from psi(R,Z) and Ra_rho_t function';
+Rminzmag=azmag;
+Rminzmag.name='Rmin(rho,t)=Rbnd_HFS from psi(R,Z) and Ra_rho_t function';
+Rmaxzmag=azmag;
+Rmaxzmag.name='a(rho,t)=Rbnd_LFS from psi(R,Z) and Ra_rho_t function';
+
+iteff=find(sspr.t>=54 & sspr.t<=69);
+% for it=1:length(sspr.t)
+for it1=1:length(iteff)
+ it=iteff(it1);
+ [r,z,psinrz]=psinrzjet(-1,sspr.t(it),3*npts,3*npts,[],[],[],0,sspr,sspi,[]);
+ ir=find(r>2.7 & r<3.3);
+ iz=find(z>-1 & z<1.5);
+ psimin=1.2*max(0,min(min(psinrz(ir,iz))));
+ [c,h]=contour(r,z,psinrz',[psimin psimin]);
+ iii=find(c(1,:)==psimin);
+ rmagpsi=mean(c(1,iii(end)+1:iii(end)+c(2,iii(end))));
+ zmagpsi=mean(c(2,iii(end)+1:iii(end)+c(2,iii(end))));
+ % make (R,Z) array on which to compute psinrz
+ nr=sspi.data(1,it);
+ rmin=sspr.data(1,it)/100;
+ rmax=sspr.data(nr,it)/100; % in [cm] but will be used in [m]
+ zeff=zmagpsi;
+ rgrid=linspace(rmin+1e-4,rmax-1e-4,3*npts);
+ zgrid=zeff*ones(size(rgrid));
+ [r2,z2,psinrz2]=psinrzjet(-1,sspr.t(it),rgrid,zgrid,[],[],[],0,sspr,sspi,[],1);
+ [a,da1]=interpos(13,r2,psinrz2);
+ ii=find(da1>0);ii=ii(1);
+ rmagmin=interp1([da1(ii(1)-1) da1(ii(1))],[r2(ii(1)-1) r2(ii(1))],0);
+ ii=find(r2>=rmagmin);
+ [RpsiLFS]=interpos(13,psinrz2(ii).^.5,r2(ii),rhopsi);
+ ii=find(r2<=rmagmin);
+ [RpsiHFS]=interpos(13,psinrz2([ii(end):-1:1]).^.5,r2([ii(end):-1:1]),rhopsi);
+ azmag.data(:,it)=0.5.*(RpsiLFS-RpsiHFS);
+ azmag.data(1,it)=0;
+ azmag.data(2,it)=max(0,azmag.data(2,it));
+ R0zmag.data(:,it)=0.5.*(RpsiLFS+RpsiHFS);
+ R0zmag.data(1,it)=rmagmin;
+ Rminzmag.data(:,it)=RpsiHFS;
+ Rmaxzmag.data(:,it)=RpsiLFS;
+ % get isoflux
+ % [c1,h1]=contour(r,z,psinrz',rhopsi.^2);
+ % for ir=1:length(rhopsi)
+ % ii=find(c1(1,:)==rhopsi(ir)^2);
+ % if isempty(ii)
+ % Riso{ir,it}=rmagmin;
+ % Ziso{ir,it}=zmagpsi;
+ % else
+ % [il ij]=max(c1(2,ii));
+ % ieff=ii(ij);
+ % Riso{ir,it}=c1(1,ieff+1:ieff+il)';
+ % Ziso{ir,it}=c1(2,ieff+1:ieff+il)';
+ % end
+ % end
+end
--
GitLab