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

+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