function [r,z,psinrz,sspr,sspi]=psinrzjet(shot,time,nrg_rg,nzg_zg,efitlab,uid,seqd,varargin);
%
% psirz : reconstruction des surfaces de flux
% 
% ce programme utilise les donnees de efit ou eftm
% 
% examples:
%  [r,z,psinrz,sspr,sspi]=psinrzjet(shot,time,nrg_rg,nzg_zg,[efitlab,uid,seq,ncont]);
%  [r,z,psinrz,sspr,sspi]=psinrzjet(50814,60,65,65,[],[],[],60,sspr,sspi); % to get plot and give sspr,sspi
%  [r,z,psinrz]=psinrzjet(50814,60.4,[3 3.2],[0 0.1],[],[],[],0,sspr,sspi,[],1); % 
%
% entrees :
% shot : numero du choc jet
% time   : time de l'analyse
% nrg_rg, nzg_zg: nb de points de la grille en r (resp. en z) sur laquelle on fait la 
%           reconstruction des surfaces de flux.
%           if nzg_zg is negative, make symmetric box around zero
%           if array, assumes rout and zout given firectly
% varargin{1}: plot option: 0: do not plot contours, >0 plot contour with varargin{1} nb of contours (60 is good)
% varargin{2}: sspr from ppf/efit/sspr, structure containing sspr.data and sspr.t
% varargin{3}: sspi from ppf/efit/sspi, structure containing sspi.data and sspi.t
% varargin{4}: deltaz
% varargin{5}: idiag: =1 : gives diag(psinrz) as output to give psi at (r,z) points
%              idiag=0 (default) full matrix as output
%
% facultatifs :
% efitlab : efit ou eftm (efitlab=efit par defaut)
% uid     : eventuellement, donnees ppf privees (uid='jetppf' par defaut).
% seq     : "sequence number" de l'uid (0 par defaut -> version la plus recente)
%
% sorties :
% r, z    : vecteurs de la grille (r,z)
% psinrz  : matrice psin(r,z) (flux normalise)
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if nargin<=4 | isempty(efitlab)
  efitlab='efit';
end
if nargin<=5 | isempty(uid)
  uid='jetppf';
end
if nargin<=6 | isempty(seqd)
  seqd='0';
end
ncont=0;
if nargin>=8 & ~isempty(varargin{1})
  ncont=varargin{1};
end

% equilibre magnetique
%---------------------
if nargin>=9 & ~isempty(varargin{2})
  sspr=varargin{2};
else
  sspr=gdat(shot,'ppf',[efitlab '/sspr?uid=' uid '+seq=' seqd]);
end
ssprs=sspr.data;
tpefit=sspr.t;
if nargin>=10 & ~isempty(varargin{3})
  sspi=varargin{3};
else
  sspi=gdat(shot,'ppf',[efitlab '/sspi?uid=' uid '+seq=' seqd]);
end
sspis=sspi.data;
tpefit=sspi.t;
if nargin>=11 & ~isempty(varargin{4})
  deltaz=varargin{4};
else
  deltaz=0;
end
if nargin>=12 & ~isempty(varargin{5})
  idiag=varargin{5};
else
  idiag=0;
end
  
[x,ind]=min(abs(time-tpefit));
sspr_t=ssprs(:,ind);
sspi_t=sspis(:,ind);

nr=sspi_t(1);
nz=sspi_t(2);
nc=sspi_t(3);
ri=sspr_t(1:nr);
zi=sspr_t((nr+1):(nr+nz));
cij=sspr_t((nr+nz+1):(nr+nz+nc));

rmin=0.01*sspr_t(1);
rmax=0.01*sspr_t(nr);
zmin=0.01*sspr_t(nr+1);
zmax=0.01*sspr_t(nr+nz);

rmin=rmin+1e-4;
rmax=rmax-1e-4;
zmin=zmin+1e-4;
zmax=zmax-1e-4;

if length(nrg_rg)==1 & length(nzg_zg)==1
  r=linspace(rmin,rmax,nrg_rg);
  z=linspace(zmin,zmax,abs(nzg_zg));
  if nzg_zg<0
    zlim=max(abs(zmin-deltaz),abs(zmax+deltaz));
    z=linspace(-zlim-deltaz,zlim-deltaz,abs(nzg_zg));
  end
else
  r=nrg_rg;
  z=nzg_zg;
end

% mapflux contruit la carte de flux psin sur (r,z)
psinrz=mapflux_fast(cij,0.01*ri,0.01*zi,r,z,idiag);

if ncont>0
  figure
  contour(r,z,psinrz',ncont);
  hold on
  [h1, h2]=contour(r,z,psinrz',[1 1],'k');
  for i=1:length(h2)
    set(h2(i),'LineWidth',2);
  end
  ss=sprintf('%.4f',time);
  title(['JET #' num2str(shot) ' t= ' ss])
  xlabel('R [m]')
  ylabel('Z [m]')
  axis equal
  axis([rmin rmax zmin zmax])
end