function [R0zmag,azmag,Rminzmag,Rmaxzmag,zmageff]=Ra_rho_t(sspr,sspi,npts);
%
% function [R0zmag,azmag,Rminzmag,Rmaxzmag,zmageff]=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, zmageff
%
% 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';
zmageff=NaN*ones(size(sspr.t));

% 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;
  zmageff(it)=zeff;
  % 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